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PREPACB 

This  dissertation  is  concerned  with  two  problems, 
relatively  independent  of  one  another.   The  first  is  the 
question  of  the  effect  of  long  range  forces  on  the  radial 
distribution  function  g  and  the  second  deals  with 
attempts  to  improve  the  integral  equations  from  which  g 
may  be  computed  for  any  potential. 

Chapter  I  presents  a  brief  review  of  the  major 
results  obtainable  from  the  cluster  diagram  expansion  of  g 
pioneered  by  Mayer  and  Montroll.   It  provides  the  background 
for  both  Parts  One  and  Two.   The  special  case  of  long  range 
forces  is  taien  up  in  the  second  chapter,  where  two  expres- 
sions are  obtained,  based  on  different  approximations, 
which  relate  a  given  "short  range"  g  and  a  long  range 
potential  to  the  total  g  arising  from  the  sum  of  the  long 
and  short  range  potentials.   These  expressions  are  then 
solved  numerically  for  the  Gaussian  model  and  comparisons 
made  to  previous  solutions.   The  solutions  are  reported  in 
Chapter  III,   Chapters  IV  and  V  are  taken  up  by  the  generali- 
zation to  two-component  systems  of  the  expressions  derived 
in  Chapter  II  and  their  numerical  solution  for  a  model  of 
an  electrolyte. 
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In  prittclple,  the  radial  distribution  function  for 
any  fluid  may  be  found  "exactly"  by  the  Monte  Carlo  (MO) 
technique,  given  the  pair  potential  of  the  fluid.   In 
practice,  however,  this  technique  is  a  voracious  consumer 
of  computer  time,  so  that  it  is  usually  used  as  a  means  of 
producing  "touchstones"  against  which  other,  more  frugal 
approximate  solutions  may  be  tested.   These  approximation 
methods  usually  result  in  integral  equations  for  g  ,  which 
have  the  added  advemtages  over  MC  that  they  provide  a  means 
of  computing  experimental  potentials  from  scattering  data 
and  may  indeed  have  analytic  solutions,  as  witnessed  by 
Werthelm's  solution  of  the  Percus-Yevick  (PY)  equation  for 
hard  spheres. 

The  diagram  development  of  Chapter  I  suggests  a 
possible  avenue  of  approach  in  the  improvement  of  the 
integral  equations  for  g  derived  there.   This  matter  is 
taken  up  again  in  Chapter  VI,  which  begins  Part  Two.   The 
upshot  of  the  arguments  presented  here  is  a  new  integral 
equation  for  g  which  contains  a  parameter,  m,  determined 
by  the  peculiar  potential  and  temperature  of  the  system 
studied.   This  equation  reduces  to  the  well-known  Percus- 
Yevick  (PY)  or  convolutlon-hypemetted  chain  (CHNC)  Integral 
equations  upon  the  proper  choice  of  the  value  of  m  .   Solu- 
tions of  the  equation  have  been  obtained  numerically  for 
the  Lennard- Jones,  Coulomb,  and  hard  sphere  potentials  and 

It 


the  results  compared  to  the  PY  and  CHNC  solutions,  as  well 
as  to  the  "standard"  solutions  available.  These  numerical 
results  constitute  Chapter  VII. 

The  numbering  of  equations  throughout  the  disserta- 
tion is  bj  chapter  and  serial  order  therein.   Thus  Bq. 
(5.11)  refers  to  the  eleventh  equation  in  Chapter  V.   The 
same  rule  applies  to  tables  and  figures. 
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PART  ONE 

PERTURBATION  CORRECTION  OF  THE  RADIAL  DISTRIBUTION 
FUNCTION  FOR  LONG  RANGE  FORCES 


CHAPTEH  I 
CLUSTBH  EXPANSION  0?   THE  RADIAL  DISTRIBUTION  FUNCTION 

Statistical  mechanics  seeks  to  deduce  the  macroscopic 
thermodynamic  behavior  of  a  physical  system  from  the  laws 
governing  the  microscopic  behavior  of  its  constituent 
particles.  This  goal  can  be  embodied  in  the  calculation  of 
a  single  quantity  from  which  all  thermodynamic  information 
can  be  obtained.  Thus,  for  a  classical  canonical  ensemble, 
this  end  is  achieved  by  a  knowledge  of  the  partition 
function  Q  defined  by 

■J    J     i=i      1=1 

where 

3  -  lAT  .  (1.2) 

Here  T  is  the  absolute  temperature  and  k  and  h  are 
Boltzmann's  and  Planck's  constants,  respectively.   The 
N-particle  Hamiltonian  H  describing  the  system  is  a  function 
of  the  3N  components  of  momentum  and  the  $N  components  of 
position. 

The  free  energy  P  of  the  system  is  then  given  by 

-PF  .  In  Q  (1.3) 
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and  all  other  thermodTnamic  quantities  can  be  obtained  from 
the  free  energy.  For  example,  the  energy  E  of  the  system 
is  given  by 

E  -  d(3P)/d0  .  (1.^) 

Another  quantity  that  may  be  used  to  compute  the 
thermodynamic  quantities  is  the  pair  distribution  function 
g(r,  tPo)  (defined  below)  if  the  particles  are  assumed  to 
interact  only  through  two-body  forces.   The  integrations 
over  the  momenta  in  Eq.  (1.1)  may  be  carried  out  immediately, 
yielding 


Q  -  Z/NI  A^^  (1.5) 


where 


and 


-3U 


3^  ^3 


e       d^r^^    ...   d^rjj  (1.6) 


A    -  h(27nnkT)"-^/^    .  (1.7) 

Here  U(r,  , . . .  ,rj,)  is  the  potential  energy  of  the  system 
and  V  the  volume.  Wo  have  assumed  that  all  particles  have 
the  same  mass  m  so  that  the  kinetic  energy  is 
T  »  5^  p.  /2m  .  Then  we  may  write 

-  0F  »  -  ln(NIA^^)  +  In  Z  (1.8) 

and,  from  Eq.  (1.4), 


4^ 

E-|NkT  +  |/.../e    U  d^r^^  . . .  d^v^   ,  (1.9) 

Implementing  the  assumption  of  padrwise  interaction, 
we  write 

U  =  L    0   (r.-,)   ,  (1.10) 

i<J    '  «J 

where 

^ij  -  l^i  -  ^jl  (^-ll) 

and  0  (r)  is  the  pair  potential.  With  this,  Eq.  (1.9) 
becomes 

E  -  I  NkT  ♦  I  ^  /  /0  ir^2^s(r^fr^)d^r^^^r2  (1.12) 

where  we  have  defined  the  pair  distribution  function 
6(rj^,r2)  by  [1] 


6(r^.r2)  -  (1  .  1)^  J. . .  J 


-3U  3       3 

e    d'^r,  . . .  d-^r^  . 
5       N 


(1.13) 

The  5-particle,  ^particle,  etc.,  distribution  functions 
can  be  similarly  defined.  In  general,  the  n-particle 
distribution  function  is  given  by 

g^^^  (r      r  )   .t    f        fe'^^   d^r       d^r 

(l.W 
where  terms  of  order  1/N  have  been  ignored. 
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?or  a  fluid,  which  shall  be  our  exclusive  concern 
here,  g  is  a  function  only  of  the  magnitude  of  the 
separation  between  particles  1  and  2;  i.e., 

eiv^tv^)    -  6(^12^   •  (1.15) 

This  has  given  rise  to  the  use  of  the  name  "radial  distri- 
bution function."  The  normalization  of  the  pair,  or  radial, 
distribution  function  is 

\     sMd^T   -  1  -  I   .  (1.16) 

Physically,  the  distribution  function  g(r)   is 
proportional  to  the  probability  of  finding  any  given  pair 
of  particles  in  the  fluid  separated  by  a  distance  r  , 
irrespective  of  the  positions  of  the  remaining  N-2  psirticles. 
It  has  the  additional  property  that  its  Fourier  transform 
may  be  directly  measured  by  x-ray  and  neutron  scattering 
experiments  [2,5]. 

Ignoring  terms  of  order  1/N  ,  g(r)   is  defined  as 

g(r)  -  ^   r...  fe-^^  d^r^  ...  d\  .      (1.1?) 

All  thermodynamic  quantities  can  be  expressed  in  terms  of 
g(r)  . 

The  computation  of  g(r)  can  be  made  more  tractable 
by  the  use  of  Mayer  f  functions,  whereby  the  right  hand  side 


of  Eq.  (1,17)  is  expanded  in  Mayer  cluster  integrals.  This 
results  in  the  well-known  expansion  of  6(r)  in  powers  of 
the  density  [^-9] : 

G(r)  «  g(r)  -  1  (1.18) 

-  f(r)  +  [1  +  f(r)]C(r)   .  (1.19) 

Here  C(r)  is  defined  by 

oe  ^ 

C(r)  -  ^   n*i;  ^^   n  f^.  d^r^  ...  d^r,^^ 
t-1  J 

(1.20) 

where 

fij  -  fir^.)   =  e  ^^  ^  1  (1.21) 

and  a  1,2  index  has  been  suppressed  on  r  .  Each  power  of 
the  number  density  n  ■  N/V  is  associated  with  several 
integrals,  denoted  by  L    .  The  integrals  may  be  represented 
by  diagrams  according  to  the  following  prescription.  For 
each  f. ,  in  the  integrand  draw  the  points  i  and  J  and 
connect  them  with  a  line.  Points  1  and  2,  which  are   not 
variables  of  integration,  are  called  the  fixed  points  and 
are  drawn  as  white  (i.e.,  unfilled)  circles.  The  t  field 
points  are  drawn  as  black  circles.  The  symmetry  number  s 
depends  on  the  number  of  field  points  (t)  and  the  diagram 
type  (a).   It  is  equal  to  the  number  of  permutations  of  the 
field  points  that  lead  to  the  same  diagram,  that  is,  to 


integrands  having  identical  products  of  fj ^   .   Equation 
(1.20)  may  thus  be  read  as  the  sum  over  all  graphs  having 
two  distinguishable  points  fixed,  the  contribution  of  each 
graph  being  divided  by  its  symmetry  number. 

The  first  few  terms  in  the  expansion  of  C(r)   in 
this  shorthand  notation  are 


C(r) 


n  CD  + 


/A*(^*/^  ^  |<Ap 


*I<J> 


>  •  •  .  (1.22) 


where  the  density  and  symmetry  factors  have  been  explicitly 
written.  Any  diagram  can  easily  be  written  out  explicitly 
in  terms  of  the  integrations  it  represents  by  applying  the 
prescription  given  above  in  reverse  order.   Thus  we  have 

<^  ■   /  ^15^5^^1^^^2^S^H  •  (1-23) 

These  diagrams  can  be  classed  into  three  groups 
according  to  their  topological  characteristics: 

(a)  A  series,  or  nodal,  diagram  contains  at  least  one 
field  point  (called  a  node)  through  which  all  paths 
connecting  1  emd  2  must  pass. 

(b)  A  parallel  diagram  contains  at  least  two  paths 
between  1  and  2  which  are  connected  only  at  1  and 
2. 


(c)  A  bridge  diagram  is  one  which  is  neither  series 
nor  parallel. 
The  sum  of  all  series,  parallel,  and  hridge  diagrams  will 
"be  denoted  by  S(r) ,  P(r),  and  B(r),  respectively.  To  the 
second  power  in  density  these  sets  are 


S(r)  »  ncr\)-^   n 
P(r)  -  I  n2  qAo 


/~\  *,^*/^ 


+  •   • 


n2  ^ 


+  o 


and 


or 


B(r)  -  5 

In  terms  of  these  sets,  we  have 
C(r)  »  S(r)  +  P(r)  +  B(r) 

G(r)  -  S(r)  +  P(r)  +  B(r) 

+  f(r)  [1  +  S(r)  +  P(r)  +  B(r)] 


g(r)e^^^^^  -  1  +  S(r)  ^   P(r)  *   B(r)   . 


+  .  .  .  (1.24) 

(1.25) 
(1.26) 

(1.27) 
(1.28) 


(1.29) 

We  will  denote  the  set  of  all  non-nodal  diagrams  in  the 
expansion  of  G(r),  Eq.  (1.28),  by  T(r) ;  i.e., 

T(r)  -  P(r)  +  B(r)  ^   f(r)Cl  +  S(r)  +  P(r)  ^   B(r)] 

(1.50) 
-  g(r)f(r)e^^^^^  *   P(r)  ^   B(r)  (1.31) 


and  thus 

G(r)  -  S(r)  +  T(r)   .  (1.32) 

It  is  important  to  note  that  the  sets  S,  P,  and  B 
do  not  tnclud©  diagrams  with  a  direct  1-2  bond.   Diagrams 
of  this  type  are  put  in  explicitly  by  the  definitions,  as 
in  the  last  term  of  Eq.  (1.28).  Thus  also  the  non-nodal  set 
T  defined  by  Bq.  (1,30)  includes  not  only  the  sets  P  and 
B  but  also  all  diagrams  with  a  direct  1-2  bond. 

Consider  a  typical  diagram  of  the  set  S  .   By 
definition  it  has  at  least  one  field  point  which  is  a  nodal 
point.  Let  the  first  nodal  point  encountered  along  a  path 
from  1  to  2  be  labelled  3.  Then  the  subdiagram  between  1 
and  3  will  have  no  nodes,  by  construction,  and  hence  will 
be  some  member  of  the  set  T  *   The  subdiagram  between  3 
and  2  may  or  may  not  have  further  nodes.   There  is  no 
restriction.  Hence  in  general  it  will  be  some  member  of 
the  set  G  .   Thus  we  can  generate  all  possible  members  of 
the  set  S  by  substituting  between  1  and  3  all  possible 
members  of  T  and  between  3  and  2  all  members  of  G  ,  That 

18 

S(r)  -  n  /  nT^^)(iir^2^d}r^    ,  (1,33) 

The  factor  n  is  introduced  because  the  point  3  is  a  field 
point  for  the  whole  diagram,  although  it  plays  the  role  of 
a  fixed  point  for  the  subdia grams.  The  same  applies  for  the 
integral  over  r,  , 


10 

with  Eq.  (1.52)  we  have  thus  an  integral  equation  for 


G(r)  , 


Q(r)  -  T(r)  +  n  /  T(r')G(|r  -  r' |  )dV  .       (1.5^) 

This  is  the  integral  equation  of  Ornstein  and  Zernioke  [10] 
which  gives  G(r)  when  the  direct  correlation  function 
T(r)  is  known.  Unfortunately,  we  can  see  from  Bq.  (1.31) 
that  T  involves  the  unknowns  P  and  B  .  We  will  see 
how  to  eliminate  the  set  P  in  the  following,  but  the 
bridge  set  will  remain  to  prevent  the  formation  of  an  exact 
integral  equation  with  only  one  unknown. 

The  Fourier  transform  of  Bq.  (1.5^)  may  be  taken, 
using  the  convolution  theorem,  with  the  result  that 

G(k)  »  T(k)/[1  -  nT(k)]  (1.35) 

where  the  tilde  denotes  the  Fourier  transform.   (See 
Appendix  A.)  Similarly,  from  Bqs.  (1.32)  and  (1.53)  we  have 
that 

S(k)  -  nT^(k)/Cl  -  nT(k)]   .  (1.36) 

Consider  now  a  typical  member  of  the  set  P.   Since 
there  is  no  integration  over  particles  1  and  2  the  diagram 
may  be  "factored"  by  cutting  at  the  two  fixed  points  and 
separating  out  the  various  independent  branches.   For 
example,  we  may  write 


d^l  [  c/N, 


11 

2 

(1.37) 


The  branch  diagrams  may  be  members  of  the  sets  S  and  B  , 
but  not  of  P  . 

We  gire  Interested  in  the  contribution  of  all  pareillel 
diagrams  with  t  branches,  which  we  will  denote  by  P^-Cr)  . 
Consider  a  typical  diagram  with  t  branches  and  let  the 
contribution  from  the  iU)  branch  be  a^(r)  ,  so  that  the 
contribution  of  the  whole  diagram  is 

Ca3_(r)][a2(r)]  -.  Ca^(r)]   . 

We  can  generate  a  subset  of  p^(r)  by  holding 
a2(r) , . . . ,a^(r)  constant  and  allowing  ou(r)   to  be  any  of 
the  diagrams  of  P(r)   and  B(r)  .   The  contribution  of 
this  subset  is  then 

[P(r)  ^   B(r)][a2(r)]  ...  [a^(r)]   . 

A  still  larger  subset  of  P^(r)  ,  which  includes  the  one 
above  as  a  subset  of  itself,  may  be  found  by  allowing  the 
second  brsinch  to  be  composed  of  any  member  of  S(r)  or 
B(r)  as  well  as  the  first,  giving  the  contribution 

I  CS(r)  +  B(r)]CS(r)  +  B(r)][a^(r)]  ...  [a^(r)]  . 

Diagrams  with  non-equal  subdiagrams  are  counted  twice  in 
the  multiplication  of  the  first  two  bracketed  terms  and  the 
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factor  1/2  removes  this  redundancy.   Diagrams  vrith  equal 
subdiagrams  are  only  counted  once  and  here  the  factor  1/2 
serves  to  account  for  the  correct  symmetry  factor.  We  may 
go  on  in  this  way,  allowing  the  first  3  branches  to  take 
on  arbitrary  members,  and  so  on,  obtaining  for  the  total 
contribution  of  diagrams  having  t  branches 


P^(P)  .  CS(r)^|  B(r)3^   ^ 


(1.58) 


For  the  parallel  set  P(r)  ,  t  may  be  any  number 
from  2  to  infinity,  so  that 


CO 


t»2 


(1.39) 


-  e^^""^  ^  ^^""^   -  S(r)  -  B(r)  -  1  .   (1.^) 
Thus,  from  Bq.  (1.29), 

g(r)eP^(^)  .  eS(^)  *  ^(r)  ^^^^^^ 


In  g(r)e 


30(r) 


»  S(r)  +  B(r)  . 


(1.^2) 


Then  substituting  Eqs.  (1.^1)  and  (l.'<-2)  into  (1.^0)  givea 
a  useful  relation  for  P  ,  viz.. 


P(r)  -  g(r)e 


3lZ)(r) 


-  1  -  In 


g(r)e 


30(r) 


(1.^3) 


which  can  be  used  to  eliminate  P(r)   from  previous  equations, 
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Thus  the  direct  correlation  function  T  "becomes,  from 
Eq.  (1.51). 

T(r)  -  G(r)  -  ln[g(r)e^^'^^^  ]  >  B(r)  .         d.'^^) 

Unfortunately  the  set  B  has  no  exploitable  properties 
which  would  effect  its  elimination  and  it  is  at  this  point 
that  approximations  must  be  introduced.   Perhaps  the  simplest 
approximation  which  can  be  made  is  to  neglect  B  altogether, 
hoping  that  it  is  small.  With  this  approximation, 

T(r)  -  G(r)  -  ln[g(r)e^^^^^  ]  .  (1.^5) 

Inserting  this  into  Eq,  (1,5^)  gives  the  convolution- 
hypemetted  chain  (CHNC)  integral  equation  [11-16] 

G(r')  -  In  g(r')e^^^^'^ 


ln[g(r)e0^(^)]  -  n 


X  G(  |r  -  r'l  )d^r'  .        (1.^) 


Another  approximation  to  B  which  has  proved  useful 
is  the  assumption  that 

B(r)  ^  -  P(r)  (1.^7) 

or 

B(r)  ♦  P(r)  «=  0  .  (1.^8) 

Neglect  of  these  two  sets  gives 

T(r)  -  g(r)f(r)e^^^''^  (1.^9) 

as  can  be  seen  easily  from  Bq.  (1.51).  Combined  with 
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Eq.  (1.5^) »  this  approximation  leads  to  the  Percus-Yevick 
(PY)  integral  equation  [17-20] 

B(r)e^^(r)  =  1  +  n  TgCr' )f (r' )e^^(^' ^ 

X  G(|r  -  r'|)d^r'  .         (1.50) 

A  third  and  older  integral  equation  for  g  which  is 
often  found  in  the  literature  is  the  equation  of  Bom-Green- 
Yvon  (BGY)  [21,22],  based  on  the  superposition  approximation 
[1].   It  has  been  found  to  be  less  reliable  than  the  PY  or 
CHNC  integral  equations  for  a  Lennard-Jones  (6-12)  potential 
by  Broyles  et  al.  [23]  and  for  hard  spheres  by  Klein  [2^,25]. 
Ve  will  not  have  occasion  to  refer  to  it  again  in  this 
dissertation. 

For  a  realistic  potential  made  up  of  a  repulsive  core 
and  an  attractive  part  such  as  the  Lennard-Jones  (6-12) 
potential,  Broyles  et  al.  [25]  have  found  that  at  relatively 
high  temperatures  the  PY  equation  is  superior  to  the  CHNC, 
while  the  results  obtained  by  Khan  [20]  suggest  that  at  low 
temperatures  the  CHNC  equation  may  be  the  better  of  the  two. 

We  will  return  to  this  temperature  effect  in  Part 
Two,  where  an  attempt  is  made  to  obtain  an  integral  equation 
capable  of  adjusting  itself  to  the  particular  potential  and 
temperature  desired. 


CHAPTBH  II 
LONG  RANGE  FORCES 

We  will  be  concerned  in  this  chapter  with  the  effect 
on  g(r)  ,  and  hence  on  the  thermodynamic  quantities,  of  a 
small  change  in  the  pair  potential,  0(r)  ,  describing  the 
system.   A  solution  to  this  problem  could  be  used  in  a 
variety  of  applications.   On  a  practical  level,  the  need 
for  a  method  to  correct  g  arises,  for  example,  in  Monte 
Carlo  calculations  of  the  radial  distribution  function, 
where  the  long  range  tail  of  a  potential  such  as  the 
Coulomb  potential  must  necessarily  be  truncated  at  some 
finite  distance.   The  effect  of  the  neglected  part  of  the 
potential  must  be  found  for  a  complete  solution  [26,2?]. 
Furthermore,  if  the  function  g  is  known  for  some 
temperature  T  ,  g  for  some  slightly  different  temperature 
t'  may  be  easily  found  by  considering  0'0  ,  p'  «  IAt' , 
to  be  a  perturbation  of  30  at  T  and  applying  the  corres- 
ponding correction.   This  obviates  the  need  of  another  long, 
iterative  solution  of  the  integral  equations  for  g  . 

From  a  theoretical  point  of  view,  the  solution  of 
this  problem  would  allow  the  effects  of  the  long  and  short 
range  parts  of  the  potential  to  be  studied  independently. 
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Thus  the  change  in  the  potential  could  be  the  addition  of  a 
wealc  long  range  interaction  to  a  short  range  repulsive  core. 
The  derivation  of  a  van  der  Vaals-like  equation  of  state  is 
based  upon  this  type  of  potential  [28,29].  A  one-dimensional 
model  of  this  type  has  been  studied  by  Kac,  Uhlenbeck,  and 
Hemmer  in  a  succession  of  papers  C 30-32]  and  was  found  to 
display  a  phase  transition  even  in  lowest  order. 

The  classic  example  of  a  long  range  potential  is 
the  Coulomb  potential  which  is,  in  principle,  of  infinite 
range.  An  explicit  approximation  to  the  radial  distribution 
function  for  a  system  of  particles  interacting  with  Coulomb 
forces  was  found  by  Debye  and  Htlckel  [55]  and  reads 


g(r)  ■  exp 
where 


pe^  -rA 


(2.1) 


\  -  [^  t;  nPe^]"^/^  (2.2) 

and  •  is  the  electron  charge. 

This  result  can  easily  be  obtained  from  the  develop- 
ment of  the  last  chapter.  Prom  Eq.  (1.^1)  we  get 

sMe^^^'^   .   eS(^)  (2.5) 

by  neglecting  B(r)  ,  as  in  the  CHNC  approximation.  The 
corresponding  direct  correlation  function  T(r)  is  then 
given  by  Bq.  (1.^5) t  which  can  be  written 
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T(r)  -  G(r)  -  ln[l  ^   G(r)]  -  00(r)  .  (2.4) 

?or  a  weak  long  range  force,  G(r)  will  be  small  over  a 
large  region  so  that  the  logarithm  in  Eq.  (2,^)  may  be 
expanded.   Thus, 

T(r)  -  -  00(r)  +  I  G^(r)  +  .  .  .  .  (2.5) 

If  only  the  first  term  is  retained  we  obtain  a  non-iterative 
solution  for  g  by  combining  Eq.  (2.5)  with 

S(k)  -  n[00(k)]2/Cl  ^  n30(k)]  ,  (2.6) 

obtained  from  Sq.  (1»56)  and  the  above  approximation  to  T  • 
Then 

S(k)  -  00(k)  -  30(k)/[l  4-  n00(k)]  (2.7) 

and 


g(r)  -  exp 


CO 


2-nrr   Jq  1  -»■  n00(lc) 


k  sin  krdk 


(2.8) 


where  the  angiilar  integrations  in  the  Fourier  transform  have 
been  carried  out.  For  the  Coulomb  potential, 

P0(k)  -  ^  n3e?^^  .  (2.9) 

(See  Appendix  A).  By  inserting  this  expression  into  Bq, 
(2.8)  and  performing  the  integration  C5^],  we  obtain  the 
Debye-Httckel  (DH)g  .  Bq.  (2.1). 
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The  two  essential  approximations  of  the  DH 
expression  are  thus 

B(r)  ^  0  (2.10) 

T(r)  ^  -  P0(r)  .  (2.11) 

Carley  [26]  has  recently  shown  that  the  DH  g  is  at  least 
as  good  as,  if  not  better  than,  g's  obtained  from  the  PY 
and  CHNC  integral  equations,  both  of  which  evaluate  many 
more  diagrams  in  the  expansion  of  g  . 

Let  us  now  consider  a  potential  0(r)  which  can  be 
written  as  the  sum  of  a  short  range  and  a  long  range  part; 
i.e., 

0(r)  -  0^^(r)  +  0^^(r)  .  (2.12) 

We  will  assume  that  g  corresponding  to  the  short  range 
potential  is  a  known  function,  which  we  call  g   .  From 
Eq.  (1.^1)  we  may  write 

g^^(r)  -  expCS^^(r)  +  B^^(r)]  ,  (2.13) 

where  the  label  "sr"  on  S  and  B  indicates  they  are  the 
sets  correspondinLg  to  0   ,  Simileirly,  for  the  complete 
potential  0  ,  we  may  write  another  equation  identical  to 
Eq.  (1.^1).   Then  dividing  this  equation  by  Eq.  (2.15)  we 
have  formally 
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g(r)  -  g'*^(r)exp[-p0'-^'(r)  +  AS(r)  +  AB(r)]     (2.1^) 
where 

AS(r)  -  S(r)  -  S^'^(r)  (2.15) 

and  AB(r)  is  similarly  defined.  Prom  Bqs.  (1.55)  and 
(1.56)  we  have  further  that 

-2  -^srp 

Ai(k)  -  '''^  W       -  ^T   W  (2.16) 

1  -  nT(k)   1  -  nT^^(k) 

1  -  n[l+nG®^(k)]AT(k) 

All  we  need  now  to  complete  the  solution  is  a  set  of 
suitable  approximations  for  AT  and  AB  .  Both  of  these 
sets  are  due  to  the  long  range  potential  since,  if  it  were 
absent,  they  would  be  identically  zero.  Hence  by  analogy 
with  the  long  range  problem  solved  by  the  DH  g  ,  we  neglect 
AB  and  approximate  AT  by  -30   .   This  leads  then  to  the 
final  expression 

g(r)  -  g^^(r)e-^^^)  (2.18) 

with 

H(k)  .  [1  ^  ^"wAi'^C^^  . 

1  *  n[l  *  nG«'(k)]00^=^(k)  ^    ^^' 
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Here  g(r)  is  given  completely  in  terms  of  the  long  range 
potential  suid  the  short  range  g  ,  which  is  the  desired  form 

C55]. 

It  is  easy  to  obtain  higher  approximations  to  AT  . 
Having  written  the  potential  as  a  sum  of  two  parts,  we  may 
write  for  the  corresponding  Mayer  f  function 

f(r)  -  f^^(r)  +  e-^^^^^'^^^f^^Cr)  (2.20) 

-  f^^(r)  +  Af(r)  .  (2.21) 

If  we  replace  each  f  "bond  in  a  typical  diagram  of  the 
expansion  of  T(r)  ,  for  example, 

o^   =  /f3_2fi5f32^^^5  '  ^^-^^^ 

S3? 

by  the  sum  of  f    and  Af  and  expand  the  products,  we 
obtain  the  following, 


^       ,•>      ,^.      A, 


(2.23) 


Here  a  wavy  line  connecting  two  points  indicates  an 

ar  -Q0^^   Ir 

f    bond  and  a  dashed  line  represents  Af  ■  e  ^   f 
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Since  we  have  assumed  that  0  ^  is  a  weak  potential  we  need 
consider  only  terms  up  to  first  order  in  Af;  that  is, 
diagrams  having  only  one  dashed  line.  Performing  the  above 
expansion  in  each  diagram  of  T  gives  rise  to  one  set  of 
diagrams  having  nothing  but  short  range  bonds,  another 
having  one  AP  bond  in  each  diagram,  and  so  on.   The  sum  of 
the  first  set  gives  immediately  T^^(r)  .  We  write 


,sr 


T(r)  -  T''-'(r)  +  0----0  +  n 


CT-  -  —to  ^   o*^'^~0 


.  A> 


+  ...   .    (2.2^) 


The  sum  of  all  diagrams  having  one  AP  bond  is  the 
desired  AT  .   The  previous  approximation  to  AT  may  be 
recovered  by  taking  only  the  single  bond  diagram  after  T 
as  AT  ,  writing  it  in  two  parts, 

0.-0  -  e-P<^''"(^)fl^(r)  (2.25) 

-  f^^(r)  +  f^^(r)f^^(r)  (2.26) 

and  further  approximating  by  neglecting  the  second  part  and 
taking 

f^^(r)  -  -  P0^''(r)   .  (2.27) 
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That  is, 

AT(r)  ^   -  P0^^(r)   .  (2.28) 

We  can,  however,  retain  the  whole  first  diagram,  giving 
instead 

AT(r)  -  e"^^  f^^(r)   .  (2.29) 

Further,  all  diagreoas  with  a  direct  Af  bond  between  1  and 
2  may  be  summed  to  give 

AT(r)  =-  Af(r)g^^(r)e^^  ^^^  (2.50) 

-  g^^(r)f^^(r)   .  (2.31) 

It  will  turn  out,  however,  that  these  higher 
approximations,  Eqs.  (2.29)  and  (2.51)i  are  not  as  satis- 
factory as  the  original  approximation,  Eq.  (2.28).  This 
will  be  discussed  further  in  the  next  chapter. 

Consider  now  the  logarithm  of  both  sides  of  Bq,, 
(2.18), 

Ing(r)  -  Ing^'^(r)  -  H(r)   .  (2.52) 

In  the  region  where  g(r)   is  close  to  one,  we  may  approxi- 
mate 

Ing(r)  ^  g(r)  -  1  (2.53) 

Ing^'^(r)  -=g^'^(r)  -  1  .  (2.3^) 


Then  Eq.  (2.52)  becomes 


,  sr 


G(r)  -  G'=*-^(r)  -  H(r) 


or 
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(2.35) 


1  -t-  nG(k) 


1  -h  iiG^^(k) 


1  +  nCl  +  nG^^(k)]3?^^(k) 


(2.36) 


G(k)  -  i 


1  -K  QG^^(k) 


T^rsv, 


:^rF, 


- 1 


1  +  nCl  -^  nG^'-(k)]00-^^(k) 


-•  (2.37) 


This  is  precisely  the  equation  obtained  by  Broyles,  Sahlin, 
and  Carley  (BSO)  by  a  method  based  on  collective  coordinates 
C36,373.  The  BSC  equation,  Eq.  (2.57),  was  found  to  give, 
in  some  cases,  anomalous  negative  g's  for  small  values  of 
r  .  The  derivation  given  above  suggests  that  it  is  not 
applicable  at  small  r  ,  since  the  substitution  of  g(r)  -  1 
for  Ing(r)  is  not  valid  in  this  region. 

If  we  now  further  restrict  the  region  of  r  under 

sr 
consideration  to  that  where  g    is  practically  one,  then 

G    in  Eq.  (2,35)  may  be  neglected,  giving 

G(r)  -  -  H(r)   .  (2.38) 

This  equation  was  obtained  by  Eemmer  as  the  first  order 
expression  for  g(r)  in  the  large  r  region  [38], 

We  note  also  that  in  the  limit  0  (r)  -»  0  bo  that 
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0"^^(r)  —  0(r)  ,  Bqs.  (2.18)  and  (2.19)  give  the  DH  g 
when  0  is  specified  as  the  Coulomb  potential. 

An  equation  for  g(r)  employing  a  PY-type  approxi- 
mation, where  AP(r)  +  AB(r)  is  neglected,  is  easily 
obtained.   From  Eq.  (1,29),  we  write 

g^^(r)e^^^''<^)  -  1  ^  S^^(r)  *  P^^(r)  ^  B^^(r)   . 

(2.39) 
A  similar  equation  holds  for  the  total  g  .   Subtracting 
one  from  the  other  then  gives 

g(r)eP<^(^)  »  g^^(r)eP^'''(^)  ^  AS(r)  >  AP(r) 

+  AB(r)    (2.^0) 

g(r)  -  g^^(r)e-P<^^''(^)  .  e'^^^^^ t^M  (2.^1) 

where  we  have  used  the  approximation 

AB(r)  =^  -  AP(r)   .  (2.42) 

The  set  AS  is  determined  in  terms  of  AT  by 
Eq.  (2,17),  as  before.   If  we  now  again  approximate  AT 
by  -  00  ^  we  have  finally 

g(r)  -  g^^(r)e-^<^^''(^>  *  e'^^^^^  [00^^(r)-H(r)]  ,(2./^3) 

where  H(r)   is  the  same  as  in  Eq.  (2.19).  Once  again 
higher  approximations  to  AT  may  be  incorporated  into  th« 
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final  result,  but  it  la  found  that  the  results  are  less 
satisfactory, 

Ve  have  found  two  expressions  which  give  g(r)  for 
the  total  potential 

0(r)  -  iZl^^(r)  +  0^^(r)  (2.4^) 

sr 
when  g    Is  known.  Using  a  CHNC-type  approximation, 

AB  =  0  ,  we  were  led  to  Eq.  (2,18),   A  PY-type  approximation, 

AB  +  AP  =*=  0  ,  gave  rise  to  Sq,  (2 A3).      In  both  of  these  H 

Is  determined  by  Eq.  (2,19). 

In  the  next  chapter  these  equations  are  tested  on  a 

model  for  which  near  exact  answers  are  available  for 

comparison. 


CHAPTBH  III 
A  mJMEHICAL  TEST:   THE  GAUSSIAN  MODEL 

Tiie  equations  derived  in  the  previous  chapter  may 
be  tested  numerically  if  an  exact  solution  of  the  system 

being  studied  is  available  for  comparison.  We  need  also 

sr 
an  exact  g    for  input,  so  that  any  deviations  of  the 

computed  g's  from  the  exact  g  will  reflect  only  the 
approximations  introduced  in  the  previous  chapter.  Both 
of  these  conditions  can  be  met  by  the  so-called  Gaussian 
model,  where  the  Mayer  f  function  is  represented  as  a  nega- 
tive Gaussian  function.   This  model  was  used  to  carry  out 
numerical  computations  of  Eqs.  (2.18)  and  (2.^3),  as  well 
as  of  the  equations  resulting  from  the  higher  approximations 
to  AT  .   The  BSO  equation  was  applied  to  the  Gaussian  model 
by  D.  D.  Carley.*  The  two  asymptotic  expressions  for  g 
obtained  by  Hemmer  [38]  have  also  been  included  in  the 
results  reported. 

The  Gaussian  model  has  been  recently  investigated  by 
Helfand  and  Kornegay  [593.   These  authors  generated  and 


The  author  is  indebted  to  Dr.  Carley  for  permission 
to  include  his  results  in  the  final  comparisons,  as  well  as 
for  providing  an  independent  check  on  much  of  the  Fortran 
programs  used  in  this  part. 
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classified  the  series  and  bridge  diagrams  having  five  or  less 

field  points  and  evaluated  these  for  the  Gaussian  model. 

This  data  wsis  then  used  to  obtain  the  coefficients  ga.(r) 

t 

in 


g(r) 


-30(r) 


oo 


.  Z  n^ 


t-1 


'g^(r) 


(3.1) 


up  to  t  »  5.   At  small  values  of  density  the  g's  obtained 
in  this  way  are  essentially  exact.  The  criterion  was  that 
n  ScC^)  should  have  an  almost  negligible  effect  on  the 
total  g(r)  . 

The  Mayer  f  function  is  given  by 

f(r)  -e-^^/*)^  (3.2) 

and  the  unit  of  distance  is  selected  so  as  to  make  the 

second  virial  coefficient  of  the  pressure  identical  with 

that  of  a  gas  of  hard  spheres  of  diameter  d;  i.e., 
^  -  d       ^oo 


/  f(r)i^dr  -  -/  r^  dr  --/  e'^^^^^   r^  dr 

Jo  Jo  Jo 


so  that 


or 


d^ 
T 


a5n^/2 


ui 


1/3 


Ti    a 


1.100  a 


(3.3) 


(3. A) 


(5.5) 
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The  reduced  units  are  then  given  bj 


r/d  ,  n  -  Nd^A  . 


(3.6) 


The  potential  corresponding  to  a  solution  g(x)  is 
given  by 

.2 


00(x)  -  -  ln(l  -  e'^*^^  ^  )   . 


(3.7) 


We  can  obtain  a  g(x)  for  a  slightly  different  potential  in 
the  following  way.  Perform  a  coordinate  transformation 


sr, 


x'  »  x/m   ,  where     0  <m^  1    .     Then  the  function     0     (x)    , 


ifSr 


30'*^(x)    -   30(x/m) 


-  -  ln(l  -  e-^-21(x/m)2^      ^ 


(3.8) 
(3.9) 


will  be  of  the  same  form  as  0(x)  ,  but  go  to  zero  sooner. 
That  is,  the  transformation  x'  »  x/m  tends  to  compress 
the  potential  in  towards  the  origin.  We  will  call  0^^(x) 
the  short  range  potential.   The  corresponding  g  (x)  is 
obtained  by  replacing  x  by  x/m  and  n  by  nm^  in 
Bq.  (3.1);  that  is, 

CO 

i.i: 


g^^(x)  .  e-P^(^/^> 


t-1 


(nm5)*g^(x/m) 


.   (3.10) 


The  long  range  potential  is  given  by  the  equation 


-Ir 


sr> 


(x)  -  30(x)  -  30*'^(x)   . 


(3cll) 


29 

The  two  quantities  on  the  right  are  computed  directly  from 
Eqs.  (5.7)  and  (5.9).   This  long  range  potential  and  the 
short  range  g  computed  from  Helfand's  data  by  means  of 
Eq.  (5.10)  are  then  used  as  input  quantities  for  the 
equations  obtained  in  the  previous  chapter,  as  well  as  the 
BSC  and  Hemmer  equations.  The  computed  g*s  can  then  be 
compared  with  the  nearly  exact  g  obtained  also  from  the 
Helfand  and  Kornegay  data  using  Bq.  (5.1). 

The  programs  were  written  in  the  Fortran  language. 
They  are  straightforward  programs  and  are  not  included  here. 

In  order  to  illustrate  the  numerical  procedure  we 

will  consider  a  particular  equation,  say  Eq.  (2.18),  with 

H(x)  determined  by  Eq.  (2.19).  For  numerical  calculation 

a  maximum  value  for  x  ,  x „  ,  must  be  chosen.   The  range 

'  max  '  ° 

between  0  and  x    is  then  divided  into  N  intervals  of 
equal  size  A  .  A  function  f(x)  is  then  considered  deter- 
mined if  it  is  known  at  the  N  +  1  points  x.  =  jA  , 

J 

J  =  0,IT  .   The  calculations  of  00(jA)  ,  00®^(jA),  00^^(JA)  ♦ 
ancL  s(jA)   are  straightforward  and  require  no  comment.   To 

compute  g  (jA)  we  first  recompute  g(jA)   for  the  density 

3 
nm   and  then  treinsform  to  a  new  coordinate  i  =  j/n  to 

obtain  g  (iA)  -  g(jA/m)  .   Depending  on  the  size  of  m  , 

i     will  reach  its  maximum  value  at  some  point  before  i 

does,  leaving  g  (lA)  undetermined  for  the  remaining  points. 

In  this  region  g  (iA)  was  assigned  the  value  one,  since 
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SI* 

in  every  case  g    had  already  settled  to  its  asymptotic 


veilue  of  one  "before  reaching  the  undetermined  region.  The 
short  range  g  was  then  evaluated  at  the  same  points  as  the 
other  functions  by  linear  interpolation  of  the  results 
computed  above. 

The  Fourier  transforms  appearing  in  the  equation  for 
H(k)  are  defined  as  follows:   For  any  function  of  the 
magnitude  of  r  ,  F(r)  ,  we  define  F(k)   such  that 

F(r) ^r  /  F(k)e^^'^  d\  (5.12) 

(2n)5J 

F(k)  -  /F(r)e"^^*^  d^r  .  (3.13) 

Carrying  out  the  angular  integrations, 

F(r)  -  ~L>  /  lEp(k)  sin  kr  dk  (3.1'*-) 

rF(r)  sin  kr  dr  .  (5.15) 

For  numerical  calculation  F(r)  and  F(k)  must  be  made 
discrete  functions.   Put 

r  -  iAr  ,  k  =  jAk  .  (3.16) 

Then  we  have 
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P(i)  -  ^2^^        L    dP(J)  sin  (ijATAk)         (5.17) 
2n  lAr  J-0 


N 

F(J)  -  ^^^7^  ^   iJ'Ci)  sin  (ijATAk)   .      (3.18) 
»0 


■^S^ioC 


It  is  shown  in  Appendix  A  that  in  order  to  conserve  the 
orthogonality  of  the  sine  functions  in  this  discrete  repre- 
sentation, and  thus  the  basic  relation  of  a  function  to  its 
Fourier  transform,  we  must  put 

ArAk  -  2k/(2N-^1)  (3.19) 

where  N  is,  as  above,  the  nvimber  of  points.  Then,  putting 

A  -  2n/(2N>l)   ,  (3.20) 

we  have  finally 

P(i)  -  — T-^^ r   Z  dF(d)  sin  (iJA)     (3.21) 

2n''i(Ar)^  j-0 


^i 


p(j)  .  i^JlUu-    )      ij.(i)  sin  (ijA)   .     (5.22) 


Eq.  (3.22)  is  used  to  compute  the  Fourier  transforms 
of  30(i)   and  G®^(i)  which  then  determine  H(J)  .   The 
discrete  version  of  H(x)  is  then  found  with  Eq.  (3.21), 
giving  the  final  answer 

S(dA)  -  g^^(jA)e-^^J^^   .  (3.23) 
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Hemmer's  small  r  solution  for  g(x)  ,  Eq.  (35)  of 
Reference  58 ♦  was  solved  in  the  following  form: 

g(x)  -  g^^(x)  -  D[g^^(x)  ^  A(x)]  (3.2^) 

where 


2it^  J    1+  nTl0-^ 


(k) 


k^dk  ,  (3.25) 


oo 


and 


7/  -  1  +  /4-n  n  /  G^^(x)  x^dx  ,  (3.26) 


.sr.    5 


A(x)  -  e-P^  <-)  Y.    ^^   -V^Cx)   .       (3.27) 
t=l 

Calculations  were  made  with  Eqa.  (2,18)  and  (2,^3), 
the  BSC  equation,  and  Hemmer's  two  asymptotic  equations,  for 
densities  of  0.1,  0.2,  0.3,  and  0.^.  With  these  densities, 
the  truncated  series  for  g(x)  yields  an  essentially  exact 
solution.  For  each  value  of  the  density,  four  values  of  the 
parameter  m  were  taken,  m  ■  0.9,  0.8,  0.7,  and  0.6, 
representing  successively  larger  perturbations  of  the 
potential,  or  more  significant  values  of  0  (x)  . 

The  results  for  n  o  0.^  and  m  =  0.9,  0.8,  0.7, 
and  0.6  sire  shown  in  Figures  3.1-3.^.  We  see  that,  as 
expected,  g  computed  from  the  BSC  equation  compares 
favorably  with  the  exact  g  in  the  large  x  limit,  but 
becomes  far  too  small  as  it  approaches  the  origin.  Hemmer's 
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two  aaymptotlc  equations  yield  good  agreement  in  the  small 
and  large  x  limits  for  which  they  are  applicable  hut  allow 
no  means  of  interpolation  for  intermediate  values  of  x  . 
Of  the  two  equations  given  at  the  end  of  Chapter  II,  that 
corresponding  to  the  PY-type  approximation  gives  a  better 
result  for  g  ,  although  an  additional  infinite  set  of 
diagrams  has  been  neglected  in  its  derivation.  As  with  the 
PY  integral  equation  itself,  this  success  is  expledned  on 
the  basis  of  cancellations  between  the  sets  P  and  B 
[20].  The  potentials  for  these  cases  are  shown  in  Figure 

5.5. 

The  results  of  other  cases  for  Bqs.  (2,18)  and  (2.^3) 
are  given  in  Figures  3.6  and  3.7  in  the  form  of  rms  deviation 
from  the  exact  g(x)  ,  where  the  rms  is  defined  by 


xrms  a 


N 


i    I      Cg^^jA)  -  6^°°^P(dA)]' 


1/2 


(3.28) 


Here  N  is  again  the  number  of  points  taken  in  the  numerical 
solution  and  A  the  interval.   The  superscripts  on  g(x) 
refer  to  the  exact  and  computed  g's.  We  note  from  Figures 
3,6  and  3.7  that  the  rms  values  obtained  for  Eqs.  (2,18) 
and  (2.^3)  are  relatively  insensitive  to  changes  in  density. 
It  was  mentioned  in  Chapter  II  that,  although  higher 
approximations  could  be  found  for  AT  ,  the  results  were 
less  satisfactory  than  those  obtained  with  AT  approximated 


3^ 

by  -30   .  The  basis  of  this  assertion  lies  in  the  compari- 
son of  results  for  the  Gaussian  model  of  the  three  approxi- 
mations for  AT  ,  Eqs.  (2.28),  (2.29),  and  (2.51).   The 
situation  is  exemplified  by  Figure  5.8,  where  the  rms 
values  for  the  equations  resulting  from  the  PY-type  approxi- 
mation and  the  above-mentioned  approximations  for  AT  are 
shown  as  functions  of  the  density  for  m  ■  0.6.   It  is  seen 
that  at  small  densities  the  higher  approximations  do  indeed 
lead  to  successively  more  accurate  g*s,  but  at  the  price  of 
far  worse  answers  at  larger  densities. 
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?ig.  5. 5. -Gaussian  model  potentials.   These  ootentials 
were  used  in  obtaining  the  g's  of  Figures  3.1-5.^, 
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Pig.  3.6. -RMS  deviations.   The  rms  deviations  from 
the  exact  Gaussian  model  g(x)  for  g  computed  from 
Eq.  (2.18)  (dashed  line)  and  Eq.  (2.4-5)  (solid  line). 
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Pig.  5. 7. -RMS  deviations.   The  rms  deviations 
from  the  exact  Gaussian  model  g(x)  for  g  computed 
from  Eq.  (2.18)  (dashed  line)  and  Eq.  (2.43)  (solid 
line) . 


4-2 


0.03r- 


rms 


0.02 


0.01 


Eq.  (2.28) 

Bq.  (2.29) 
Sq.  (2.31) 


Fig.  3. 8. -Density  dependence  of  the  RMS 
deviations  from  the  exact  Gaussian  model  g(x)  for 
g's  computed  by  neglecting  AP  +  AB  and  approximating 
AT  by  Eqs.  (2.28),  (2.29),  and  (2.31).   The  value  of 
a  is  0.6. 


CHAPTER  IV 
GENERALIZATION  TO  A  TWO-COMPONENT  FLUID 

The  equations  of  Chapter  II  may  be  readily 
generalized  to  a  two-component  fluid,  such  as  an  electro- 
lytic solution  or  a  fused  salt.  For  such  a  system,  three 
radial  distribution  functions  axe  in  general  necessary. 
If  the  fluid  components  are  labelled  1  and  2,  the  required 

distributions  are  Bii%   S22'  ^^^   ^12  °  ^21* 

Ve  will  assume,  as  in  Chapter  I,  that  the  potential 
U  of  the  system  may  be  written  as  a  sum  of  pair  inter- 
actions, as  follows: 

^^^1 V    -         Z.       %C^ij)    ^  L      <^22^^iJ^ 

l£i<d£Nj_  N^+l£i<J£N 

*    L    ^i2^^ij)    .  c^-^) 

J-N^+1,N 

Here  0,,  is  an  interaction  between  molecules  of  component 
1,  022  hetween  those  of  component  2,  and  0,  p  a  "cross" 
interaction  between  members  of  components  1  and  2.   The 
numbering  of  component  1  runs  from  1  to  N,  and  of  component 

*^3 


i^4 


2  from  N,  +  1  to  N,  the  total  number,  where 


N  -  Nj_  +  N2  .  ('<-.2) 

Again  we  assume  that  we  may  write 

where  J25^g  is  the  short  range  interaction  between  components 
a  and  &   ,  and  so  on.  In  terms  of  diagrams,  the  pair 
correlation  fxmction  G-(r)  is  again  described  by  an  equation 
like  Bq.  (1.28),  where  now  three  different  G's  must  be 
distinguished.  Ve  write 


»a3<^>  •  ^aP^^)  *  ^^  *  ^aP^^^^^^aP^^)  ^  ^aP^^) 


+  B^p(r)]  ,  a,P  -  1,2         {^A) 


where  S  n  is  the  set  of  all  nodal  diagrams  which  begin 
at  a  fixed  point  belonging  to  component  a  and  end  on  a 
fixed  point  belonging  to  component  P  .   Similar  definitiona 
hold  for  the  other  sets.  Note  that  these  diagrams  are  not 
the  same  as  the  diagrams  of  Chapter  I.  While  the  topological 
shapes  are   the  same  for  both,  the  addition  of  a  second 
component  to  the  fluid  allows  many  distinct  diagrams  to  be 
made  up  from  one  diagram  of  the  single  component  type. 

In  the  following  a  and  P  will  always  be  under- 
stood to  take  on  the  values  1  and  2. 
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Direct  correlation  functions  '^ncS^-'^^   ^^®  defined 
in  complete  amalogy  with  the  one-component  fluid,  viz., 

but  now  the  expressions  for  the  series  diagrams  become 

^c^&M   -      ^  n^  /  T^(r')G^3(|r-?'|)dV  .     (4.6) 
\-l,2  J 

Here  n,   and  np  are  the  number  densities  of  components 
1  and  2,  respectively.   Equation  (4.6)  expresses  the  fact 
that  though  the  two  fixed  points  of  the  set  S  g  are 
necessarily  members  of  the  a  and   3  components,  one  from 
each,  the  intervening  field  points  may  be  of  either 
component  1  or  2,   The  field  point  r'  in  Bq.  (4,6)  is 
occupied  by  a  member  of  component  \  .   It  is  treated  as  a 
fixed  point  by  the  diagrams  in  the  sets  T  ,   and  G,  « 
and  the  fact  that  it  is  a  field  point  of  the  whole  diagram 
accounts  for  the  presence  of  the  n,   as  well  as  the 
integration  over  r  '  . 

The  differences  in  the  diagrams  may  be  simply 
illustrated  by  considering  the  coefficient  of  the  first 
power  of  n  in,  say,  gTi(r)  .  Here  we  have 


o^'^o  +    no  ©■^''^o 


n-i  \j       \j     T         ixp 
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where  the  circles  denote  component  1  and  the  squares  com- 
ponent 2.   In  the  one  component  fluid  there  was  only  one 
diagram  for  this  coefficient.   [See  Bq.  (1.22).] 
Equation  (^.^4-)  may  be  rewritten  to  read 

30^B(r) 

The  parallel  set  may  be  eliminated  as  before  to  yield 

ga3(r)e   ^^    -  e  ""^       ^*^     .  (4.8) 

That  this   is  so  can  be  made  plausible  by  considering  the 
first  few  terms  of   a  specific  g,    say  g-iT  .      Then  we  have 

n,    o-'^'^  >o  +  n«  cf\      +  n,       cr  >o 


S^-j^(r)    »  n^  o^      >o+  np 


1 


12  "*"     12  x>  +  np    o^  td 


+  n,      o^^^^'^o  +     ^1^2  ^^^^^^  +  ^1^2    ^^^^^^^ 


+  n2      o^^^'>o+        Hi      cr^^^D  +  n,  np  o"''^^"^^ 

+  nj^n2   0^'*''^+     n2^c/*^^+    .    .    .  (4.9) 


and 


^7 


(4.10) 

where  the  circles  and  squares  again  refer  to  components  1 

and  2,  respectively.   [Of.  Bqs.  (1.24)  and  (1.26)].  Then 

the  equation 

oo 

Pa0<-)  ■  Z  irrfSapC^)  *   B^g(r)]»  (4.11) 

m«2 

[see  Bq.  (1.59)]  becomes,  for  ^^^^^(r)  , 
P^3^(r)  -  I  ^x^M   +  S^^(r)B^3_(r) 

+  \  \^W   +  ...  (^.12) 

■  I  ^1^  <^  -^  ^1^  <^  *  I  ^2^  <^  *  •  •  • 

(^.13) 
from  Eq.    (4.9).      These  terms  satisfy  the   definition  of  the 
parallel  diagrams  for  the  two-component   system. 

For  the  short  range  potentials  alone,    an  equation 
identical   to  Bq.    (4.8)  may  be  written,   %ri.th  all  quantities 
labelled  with  an  "sr"   superscript.      Then  the  ratio  of  the 
two  equations  gives,    as   in  Chapter  I, 

(4.14) 


where 


and  a  similar  definition  holds  for  AB  «  . 

ocp 

Before  proceeding  to  the  evaluation  of  the  quanti- 
ties AS  Q     and  AB  q  ,  we  need  to  obtain  explicitly  the 

ap  ap 

relations  "between  the  Fourier  transforms   of  the  three     G  « 

ap 

and  the  three  T  g  implied  by  Eqs.  (^.5)  and  (^.6).  By 
taking  the  Fourier  transforms  of  these  equations,  with  the 
aid  of  the  convolution  theorem,  we  get  the  following: 

\-l,2 


These  equations  may  be  solved  simultaneously  for  either 
they  give 


Go  in  terms  of  T  □  ,  or  vice  versa.   In  the  first  case 
ap  ap 


[l-npTp2(lc)]T.,(k)  +  npT,p^(k) 

G-,(k)  -  :r^-^ ^: ^  ^^  ^  M C^.l?) 

^^      [l-n3^T^^(k)][l-n2T22(k)3  -  n3^n2T^2^^) 


[l-n,T,,(k)]Tpp(k)  +  n,T,p^(k) 

Qpp(k)  -  =-i-ii ^1 ^  ^^     ^     '^ (^.18) 

^"^      [l-n3_T^3^(k)][l-n2T22(k)]  -  n^n2T^2^^^ 


T,p(k) 

0,p(k)  »  ^ ^^-:z :r— 2 .(^.19) 

^'^  [l-n3^T^3^(k)][l-n2T22(k)]  -  n3^n2T^2^^) 
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For  the  inverse  solution  we  obtain 

[l+npGpp(k)]G,,(k)  -  npG,p^(k) 

T-,(k)  -  zr-^-^^ ~ ^       ^     '3 ^^.20) 

^^  [l+n3^G3^j^(k)][l+n2G22(k)]  -  n3^n2G3^2^^) 

[l+n,G,,(k)]Gpp(k)  -  n,G,p^(k) 

J   (k)  -  ^  ^= -    (^,21) 

22      [l+n;|_G3_3_(k)][l+n2G22(lc)]  -  n^n2G^2  (^^ 


G,p(k) 

rp  (y\    ,  i£ .(^1.22) 

^2      [l^.n3_G^^(k)][l+n2G22(k)]  -  ^^^2^12^(1^)' 


Prom  the  above  it  is  easy  to  show  that  the  following 
are  true: 

1  +  n^^G^^^Ck)  .  Cl-n2T22(k)]  /  D^(k)  (^.23) 

1  +  n2G22(k)  -  [l-n^T3^3^(k)]/D^(k)  (4.2^) 

1  -  n^^T^^Ck)  .  [l^.n2G22(k)]/Dg(k)  (^.25) 

1  -  n2T22(k)  -  [l+n^G^^^Ck)  ]  /  Dg(k)  (^.26) 

Dg(k)  -  l/D^(k)  (^.27) 


where 

Dg(k)  -  [l+n^G]^^(k)][l+n2G22(k)]  -  n3^n2G\2^(l'^)   (^.28) 


7;     2, 


D^(k)  -  [l-n3_T3_3^(k)][l-n2T22(^)^  -  n^n2T3_2^(k)  .(^.29) 
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We  now  return  to  the  problem  of  finding  expressions 
for  the  AS^g  similar  to  Bq.  (2.17).  We  will  do  the  case 
of  ^^11    •   The  others  are  similar. 

By  definition,  we  have 

Aa^l(k)  -  G^3_(k)  -  \^W   -  mW   +  f^i(k)     (^.50) 

-  \^W   -  mW   -  ^\^W    .  (^.31) 

Prom  Eq.  (^.25),  we  write 

1  +  ih.^ii(k)  »  [l-n2T|^(k)  -  n2AT22(k)]  /D^(k)   (^.32) 
where  0+;^^'^  ^^   written  in  an  expanded  form  by  replacing 

D^(k)  -  Cl-n3^T^J^(k)][l-.n2T|^(k)]  -  n^n2T^|^(k) 

-  n3^AT^3^(k)[l-n2T||(k)]  -  n2AT22(k)  [l-n^TjJ(k)] 

-  2n^n2T^2^k)AT^2^^^  '  nj^a2[AT3^3^(k)AT22(k) 

-  ATj^g^Ck)]  .  (^.33) 

We  divide  top  aind  bottom  of  the  right  hand  side  of  Eq. 
(^.52)  by  D^^(k)  ,  where  all  the  sets  of  Eq.  (^.29)  acquire 
an  "sr"  label,  and  use  Eqs.  (^.23),  (4-. 2^),  and  (^.2?)  to 
obtain 


51 

1   +  n^G^^^Ck)   -    [l+n^^G^JCk)   -  n^^T^^WD^'^W}  /  C(k) 

where 

C(k)    -  1   -  n^AT^3^(k)[l+n3^Gj[(k)]    -  n2AT22(lc)  C1+ii2G|2(1c)] 


-  2nj^H2^12^^^^^12^^^   *  n^n2[AT^j^(k)AT22(lc)   - 

Putting  Gj^.j^(k)   from  Bq.  (^^.3^)  into  Eq.  (^.51)  gives, 
finally,  after  some  straightforward  algebra, 

£^S^^W    -   [l+n^G^^(k)]^AT3^l(k)  +  n2^G^2^(k)AT22(k) 

+  2ii2G^2^k)AT3^2^^^'^^-'^^l^^^^ 

-  n^n2^^^11^^)^^22(k)-AT^2^(k)][l4.ii^G^^(k)] 


X  D^^(k) 


/  C(k)  .     (^.56) 


Note  that  if  1I2  »  0  or  AT-,2  ■  ^^^22  *  ^  ♦  *^®  above 
reduces  to  Eq.  (2.17). 
We  will  take 


AT^pCk)  .  -  B0llW  (^.37) 

as  in  Chapter  II.  For  purposes  of  simplification,  we  will 
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restrict  the  following  dtacussion  to  those  long  range 
potentials  that  satisfy  the  relation 

AT^3^(k)AT22(k)  -  AT^2^(k)  -  0   ,  (4.58) 

Since  most  of  the  two-component  fluids  which  would  be  of 
interest  here  are  made  up  of  charged  particles,  such  as 
fused  salts  and  electrolytic  solutions,  we  may  further 
restrict  the  equations  above  by  considering  only  those 
long  range  potentials  that  satisfy  the  relations 

0llW   -  0llM   -  -  0llM   -  0^^(r)  (4.39) 

which  hold  for  the  Coulomb  potential.   Then,  with  Eqs. 
(4.57)  and  (4.59),  Eq.  (4.36)  becomes 

[lH-w,(k)]23?^^(k) 

ASTT(k) ::=-i :=T-- ^ :^y- 

^^       1  +  n^[l+Wj^(k)]3r^(k)  +  n2Cl+W2(k)]30-^^(k) 

-  a\iW  C^-.^o) 

where 

w^(k)  -  n3^GjJ(k)  -  n2Gj2^k)  (4.41) 

W2(k)  -  n2G|2(k)  -  nj^G^2(^)   •  (4.42) 

We  have  finally,  then,  by  neglecting  AB^n(r)   in 
Eq.  (4.14)  , 
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fir     -^ll^^^ 

gll(r)  -  sllMe     ^^  (^.^3) 

where 

[l+w.(k)]230^^(k) 
In  a  completel7  similar  manner,  we  obtain 

gl2(^)  -  812^^)®  (^•''^^ 

where  Hp2  is  obtained  by  interchanging  1  and  2  in  Eq. 
(^.^6),  and  H,2  ^s  gotten  from  the  following  equation, 

[l4ir  (k)  ]  [l^-Wp(k)  ]  0-^w 

H,  o(k)  -  ^ ;=r=^ = ^^rz: • 

^2      1  +  n^[l+vij_(k)]30^^(k)  +  n2[l4W2(k)]P0^^(k) 

It  is  important  to  note  that  the  above  equations  for 
H  g  apply  with  the  restriction  on  the  long  range  potentials 
imposed  by  Eq.  (^.59).   The  general  H^g  may  be  obtained 
from  Eq.  (^.36)  and  the  two  other  equations  similar  to  it 
for  AS22  and  AS-,  2  • 

The  Pl-type  approximation  may  be  applied  as  in 
Chapter  II ,  giving 
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from  Eq.  (4.7)  and  its  short  range  equivalent.  Applying 
the  same  restrictions  on  the  long  range  potential  as  above, 
we  have  the  final  expressions  for  the  g's  , 

sr    '^^^""M         -30(r)   ,_ 
gll(r)  -  gii(r)e        +  e      C30^''(r)  -  E-^^M:\    , 

(4.49) 

sr,  s  -0^^''(^)    -30(^)   ir 
g22(r)  -  g|2(r)e        +  e      [30^''(r)  -  H22(r)]  , 

(^.50) 
and 

(^.51) 

where  H  q  is  defined  as  above. 

cxp 


CHAPTER  V 
A  NUMERICAL  APPLICATION:   ELBCTROLTTIC  SOLUTIONS 


The   distribution  functions  g 


g 


.  and  g, 


describing  an  electrolytic  solution  may  be  found  using  the 
results  of  the  previous  chapter.  We  will  take  as  a  model 
of  the  electrolyte  a  system  of  nonpolarizable  ions  of  charge 
+  ze  interacting  with  a  hard  core  short  range  and  a  Coulomb 
long  range  potential.  Let  R   and  R   be  the  hard  core 
diameters  of  the  positive  and  negative  ions,  respectively. 
The  interaction  potentials  are  then  the  following: 


(5.1) 
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0^^(r)  -  z^e^/Dr 


0lf(r)  -  -  0^f(r)  , 


(5.5) 
(5.6) 


where  D  is  the  dielectric  constant  of  the  solvent.  The 
last  two  equations  show  that  this  model  satisfies  Bq.  (^.^1), 
so  that  we  may  use  the  resulting  simplified  equations  of 
the  last  chapter. 

The  short  range  distribution  functions  maj  be  most 
easily  obtained  from  the  Percus-Tevick  integral  equation. 
Lebowitz  [^0]  has  recently  solved  this  equation  analytically 
for  binary  mixtures  of  hard  spheres,  generalizing  a  method 
due  to  Wertheim  [^1,^2],   In  Lebowitz' s  solution,  the 
direct  correlation  functions  of  the  hard  sphere  mixture 
are  given  by  the  following  expressions: 

(5.7) 


T  (r)  -  < 


-[a  -♦■  b  r  ■♦-  dr-^]  ,  r  ^  R 
±    ±  ± 

0  ,  r  >  R^ 


T_^_(r) 


-a 


«a^  -  —[b+^Xdx+dx^^]  »  X  ^  r  £  R^_   (5.8) 
0  ,  r  >  R^^ 


where 

X  -  r  -  \  . 

The  constants  a,b,a,b__,b,d,  and  \  are  defined  in 
Appendix  B.   We  have  assumed,  for  def initeness,  that 
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R  ^  H   in  vn^itiag  the  equations  above.   The  equations  for 
R  >  R_  can  be  easily  found  by  interchanging  all  +  and  - 
labels. 

With  the  direct  correlation  functions  on  hand,  th« 
corresponding  radial  distribution  functions  are  found  from 
Eqs.  (^.19)-(^.21),  which  relate  the  Fourier  transforms. 
Then  the  total  g's  may  be  computed  by  applying  the  results 
of  the  last  chapter,  Eqs.  (^.^5),  (^.^7),  (^.A-8),  and 
(^.51)-(^.53)»  vd-th  the  H's  determined  from  Bqs.  (^.46)  and 
(^.^9). 

This  procedure  gives  a  non-iterative  scheme  for 
computing  the  radial  distribution  functions  describing  an 
electrolytic  solution,  in  so  far  as  the  model  chosen  repre- 
sents such  a  physical  system,   A  Monte  Carlo  solution  of 
this  model  has  recently  been  carried  out  by  Shaw  [^^-3].  We 
have  computed,  for  purposes  of  comparison,  the  ionic  radial 
distribution  functions  for  the  same  input  conditions  as  the 
Monte  Carlo  solution,  although  it  is  pointed  out  by  Shaw 
in  his  concluslonfl  that  the  ensemble  of  100,000  configu- 
rations generated  for  the  Monte  Carlo  computation  "yields 
radial  distribution  functions  still  showing  fluctuations 
large  enough  to  preclude  strictly  quantitative  comparison 
with  theory."   That  is,  the  Monte  Carlo  solution  will  not 
afford  more  than  a  good  qualitative  test,  due  to  insufficient 
computer  runs. 
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Shaw  considered  an  electrically  neutral  collection 

o 

of  cations  and  anions  of  equal  size  (diameter  ^.5  A)  and 
opposite  charge  (>  one  electronic  charge)  in  a  dielectric 
continuum  with  the  pemnittivity  of  water  at  25"  C.   The 
density  was  that  of  a  0.00911  molar  electrolyte,  that  is, 

6  0  2 

a  number  density  of  2.7^  x  10   ions/A^  for  each  ionic 
species. 

These  parameters  contain,  from  our  point  of  view, 
two  unfortunate  choices,  tending  to  limit  the  model  features 
which  may  be  tested.   First,  the  distribution  functions 
ft         and  g    will  be  identical,  due  to  the  assumption  of 
equal  hard  core  diameters  for  the  ions.   The  Monte  Carlo 
computation  thus  gives  no  idea  of  the  effects  of  different 
hard  core  ion  sizes.   Second,  the  very  low  ion  concentra- 
tions puts  the  problem  in  a  region  where  practically  any 
theory  will  tend  to  give  qualitatively  correct  solutions. 
That  is,  the  approximations  usually  introduced  to  obtain 
tractable  equations  all  tend  to  be  better  at  small  densities. 
Solutions  at  higher  densities  would  "strain"  these  approxi- 
mations more,  giving  a  better  idea  of  their  region  and 
degree  of  usefulness.   At  the  present  concentration,  the 
effect  of  the  hard  cores  is  limited  to  providing  a  cut-off 
for  the  g's  in  the  small  r  region,  the  short  range  g's 
being  nearly  one  everywhere  outside  the  haird  core  diameters. 
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The  results  of  the  calculation  are  shovm  in  Figure 
5.1.   The  curves  labelled  DH  were  computed  by  Shaw  [^33 
using  a  somewhat  modified  Debye-Htickel  expression  obtained 
by  him,  where  an  effort  was  made  to  avoid  the  usual  DH 
assumptions,  viz.,  (1)  that  only  the  interaction  between 
oppositely  charged  particles  is  important,  and  (2)  that  the 
total  number  of  ions  of  either  charge  type  is  very  large. 
The  comparison  with  the  Monte  Carlo  solution  is,  as 
expected,  not  conclusive,  although  it  can  be  said  that  the 
present  solution  is  at  least  qualitatively  correct.   The 
g's  computed  by  using  a  "Pl-type"  and  a  "CHNC-type" 
approximation,  as  discussed  previously,  both  gave  essentially 
the  same  answers  and  a  single  curve  has  been  drawn  to 
represent  both. 
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PART  TWO 

inPROVMENT  OF  THE  INTEGRAL  EQUATIONS  FOR  THE  RADIAL 
DISTRIBUTION  FUNCTION 


CHAPTER  VI 
A  NEW  INTEGRAL  EQUATION 

Let  us  recapitulate  briefly  the  results  of  the 
diagrammatic  analysis  of  Chapter  I.   We  found  there  that 
the  radial  distribution  function,  g,  for  a  fluid  composed 
of  molecules  interacting  through  two-body  forces  alone, 
may  be  obtained  with  the  aid  of  the  Ornstein-Zemicke 
integral  equation,  Eq.  (1.5A-),  when  the  direct  correlation 
function  T  is  known.   In  addition,  T  could  be  written  as 
in  Eq.  (1.51),  with  two  excess  unknowns,  or  as  in  Eq,  (1.44), 
with  only  one.  This  last  equation  was  the  limit  of  the 
ability  of  the  diagram  technique,  and  a  closed  set  of 
equations  for  g  was  not  achieved. 

This  impasse  has  necessitated  the  introduction  of 
approximate  integral  equations,  foremost  among  which  are 
the  PT  and  the  CHNC  equations.  We  saw  that  they  could  be 
viewed  as  two  different  approximations  to  the  set  B  , 
namely 

B(r)  -  -  P(r)     (PY)  (6.1) 

and 

B(r)  -  0         (CHNC)   ,  (6.2) 

where  P  is  the  set  of  parallel  diagrams, 
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Numerical  solutions  of  the  PT  and  CHNC  equations 
have  shown  that 

1)  the  PY  equation  Is  superior  to  the  CHNC  equation 
for  the  hard  sphere  potential  [253 

2)  the  PT  equation  Is  superior  to  the  CHNC  equation 
for  the  Lennard-Jones  (6-12)  potential  at  rela- 
tively high  temperatures  [25] 

5)  the  CHNC  equation  is  probably  superior  to  the  PY 
equation  for  the  Lennard-Jones  potential  at 
relatively  low  temperatures  [20] 

^)  the  CHNC  equation  is  superior  to  the  PY  equation 
for  the  Gaussian  model  [^^]. 

Khan  [20]  has  shown  how  these  results  may  be  quali- 
tatively understood  by  examination  of  the  first  terms  In 
the  diagram  expansion  of  the  sets  P  and  B  ,  as  done 
below.   [See  eqs.  (1.25)  and  (1.26)].  We  shall  call  these 
P2  and  Bp ;  1 . e . , 

At  sufficiently  low  densities,  Pp  and  Bp  will  describe 
the  approximate  behavior  of  their  respective  total  sets. 
We  will  suppose  in  what  follows  that  this  is  the  case. 
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The  difficulty  with  the  set  B  is  already  evident 
in  Bp  .   The  presence  of  the  cross  bond  connecting  the 
two  field  points  prevents  the  factorization  of  the  diagram 
in  either  direct  space,  as  with  the  P  set,  or  in  Fourier 
transform  space,  as  with  the  series  set  S  .   The  PT  and 
CJHNC  equations  circumvent  this  difficulty  by,  in  effect, 
substituting  a  constant  for  the  cross  bond  in  Bp  >  equal 
to  -1  and  0  ,  respectively.   Consider  now  the  forms  of 
the  Mayer  f  bonds  for  the  various  cases  mentioned  above. 
These  are  shown  in  Figure  6.1. 

The  f  bond  for  the  haird  sphere  potential  is  precisely 
-1  up  to  the  hard  sphere  diameter,  as  approximated  by  the  PT 
equation.  Beyond  this  point  it  is  zero,  but  the  PT  approxi- 
mation is  not  ruined,  because  the  presence  of  other  bonds 
connecting  the  two  field  points  in  B2  tends  to  cut  off 
the  integrals  anyway,  taking  over  the  task  of  the  real 
f  bond.   Clearly  for  this  case  the  CHNO  approximation  is  a 
poorer  one,  since  it  approximates  the  cross  bond  by  zero 
everywhere . 

The  f  bond  for  the  LJ  high  temperature  case  is 
similar  to  the  hard  sphere  f  ,  the  effect  of  the  attractive 
well  of  the  LJ  potential  being  minimized  by  the  high 
thermal  energy  of  the  particles  of  the  system.   The  same 
arguments  apply  as  did  for  the  hard  sphere  case,  and  the  PT 
approximation  would  appear  to  be  more  reasonable. 
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The  LJ  low  temperature  case  is  characterized  by 
the  large  peak  In  the  f  bond,  resulting  from  the  potential 
well.   In  this  case,  as  the  diagram  Bp  is  integrated  over 
the  separation  distance  between  the  field  points,  the 
negative  contribution  picked  up  at  small  r  will  tend  to 
be  cancelled  by  the  positive  peak  in  f  at  larger  r  , 
resulting  in  a  smaller  net  contribution  than  for  the  high 
temperature  case.   The  net  effect  will  depend  on  the 
particular  temperature,  but  as  the  temperature  is  lowered 
the  change  will  be  toward  the  CHNC  approximation,  as  noted 
by  Khan  [20]. 

That  the  effect  of  the  cross  bond  is  far  from  being 
-1  everywhere  is  clear  in  the  case  of  the  Gaussian  f  ,  and 
here  the  CHNC  approximation  is  the  more  reasonable  of  the 
two. 

The  above  analysis  is,  of  course,  simply  qualitative, 
and  it  can  say  nothing  of  the  expected  behavior  at  high 
densities.  Nonetheless,  the  following  procedure  suggests 
itself.   Approximate  the  set  B  by 

B(r)  -  mP(r)  (6.5) 

where  m  is  a  constant  and  P  is  the  parallel  set,  and  let 
m  be  chosen  so  as  to  reflect  the  peculiar  circumstances  and 
potential  of  the  system  being  studied.   For  the  most  part, 
m  will  be  expected  to  be  between  the  CHNC  and  PY  values  of 
0  and  -1  ,  respectively. 
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Tlie  problem  is,  of  course,  to  find  a  consistent 
method  of  calculating  m  for  any  given  potential  and 
thermodynamic  conditions.   Several  possibilities  may  be 
considered. 

An  appealing  and  quite  general  approach  would  consist 
of  exploiting  the  existence  of  the  two  expressions  for  the 
pressure.   These  are,  from  Appendix  C, 

^   »  1  -  t  f r30'(r)g(r)d5r  (6.6) 


and 


-  V  ll   »  [1  ^  n  /  G(r)(^]  /nkT  .  (6.7) 

It  is  shown  in  Appendix  C  that  these  two  equations,  with 
the  approximation  of  Eq.  (6.5),  give  the  following  condition 
for  m  : 

m  -  -  1  -  C3_/C2  (6.8) 

where 

C^^  -  /[f(r)S(r)  +  I  rf'(r)S(r)  *  |  nrf '(r)^|i^]d5r 

(6.9) 
and 

C2  »  /  [P(r)  +  f(r)P(r)  ^   |  rf '(r)P(r)  ^ 

1  arf'(r)^2i£l]d5r    (6.10) 
b        fjn       • 
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The  sets  S(r)   and  P(r)   in  the  definitions  of  C,   and 
Gp     may  then  be  obtained  approximately  from,  say,  the  PT 
integral  equation. 

The  above  procedure,  however,  has  serious  computa- 
tional drawbacks.   In  the  first  place,  the  derivatives  with 
respect  to  the  density  are  not  readily  known  and  it  has 
proved  difficult,  at  least  for  the  LJ  potential,  to 
obtain  this  quantity.   Secondly,  to  obtain  m  by  the  above 
procedure  would  entail  solving  another  integral  equation, 
such  as  the  PT  equation,  as  a  first  approximation.   Thus 
the  advantage  of  generality  entails  a  corresponding  compu- 
tational disadvantage. 

A  more  limited  condition  for  m  is  derived  from  the 
virial  expansion  for  the  pressure  and  has  been  used  by 
Hutchinson  and  Rushbrooke  [^5]  in  computing  virial  coeffi- 
cients of  a  hard  sphere  gas.  We  consider  the  virial  diagram 
(where  all  variables  are  integrated  out)  derived  from  Bp 
and  replace  the  cross  bond  by  a  constant;  i.e.,  we  write 


<I>  •"  <> 


(6.11) 


The  constant  m  then  corresponds,  in  some  sense,  to  an 
average  value  of  the  cross  bond.   It  is  a  parameter  dependent 
on  the  particular  potential  and  temperature  used  to  evaluate 
the  virial  diagrams  of  Eq.  (6.11).   Writing  this  equation 
out  explicitly,  we  have 
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/ 


^13^52^1^^2^5^'^^^1^V^^3'^^^^ 


m  »  — r  I- w. ,         (6.12) 


T 


■13  52"  14" ^2^  "1^  "2"  "3  ^ 


This  is  the  procedure  that  has  been  actually  followed, 
with  the  results  given  in  the  next  chapter.  It  has  the 
advantage  of  simplicity.   Once  a  potential  is  chosen,  m 
may  be  immediately  computed.   It  is  the  same  for  any  given 
isotherm  and  is  independent  of  density,  so  that  relatively 
few  calculations  of  m  will  be  needed  for  einy  given  problem. 

Representative  values  of  m  ,  corresponding  to  the  f 
functions  illustrated  in  Figure  6.1  as  well  as  the  Coulomb 
potential,  are  given  in  Table  6.1.   They  are  also  shown  in 
Figure  6.1  for  the  f  functions  given  there.   These  values 
of  m  tend  to  confirm  the  discussion  made  above.  We  see 
that  for  the  hard  sphere  and  LJ  high  temperature  cases, 
m  is  closer  to  the  PT  value  of  -1  than  to  the  CHNC  value  of 
0  .   On  the  other  hand,  the  m's  for  the  LJ  low  temperature 
case  and  the  Gaussian  model  favor  the  CHNC  approximation. 
[The  Coulomb  m  was  computed  here  using  an  exponentially 
shielded,  or  Debye-Htickel ,  potential  in  order  to  avoid 
divergencies.]  The  effect  of  temperature  on  m  for  the 
LJ  potential  is  shown  in  Table  6.2. 

The  new  integral  equation  to  be  solved  for  g  is 
then 
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g(r)«^^^^^  -  1  -  (l4.in)P(r) 

^  n  /  [g(r')f(r')e^^^'''^+  (l4-m)P(r ')  ]G(tr-r' |  )d5r'  , 

(6.13) 
where  P  is  given  in  terms  of  g  by  Eq.  (1.^3)  and  m 
is  defined  bj  Bq.  (6.12). 

In  the  next  chapter  we  discuss  the  numerical 
techniques  used  in  solving  Eq.  (6.15)  and  display  the 
computed  results  for  the  hard  sphere,  LJ  ,  and  Coulomb 
potentials. 
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f(r) 


LJ  (Low  Temp.) 
Gaussian  Model 
LJ  (High  Temp.) 
Hard  Sphere 


-1.0 


Fig.  6.1. -Mayer  f  bonds  for  various  potentials. 


71 


TABLE  6.1 
REPRESENTATIVE  VALUES  OF  THE  PARAMETER  m 


Hard 
Spheres 

High 
Temp. 

Low 
Temp. 

Coulomb 

Gaussian 

-0.73* 

-0.81^ 

-0.^5 

-0.37 

-0.55 

*From  Hutchinson  and  Rushbrooke,  Reference  '4-5. 
Prom  Barker  and  Monaghan,  Reference  ^6. 


TABLE  6.2 

EPPBCT  OP  TEMPERATURE  ON  m  POR  THE  LJ  POTENTIAL 

kH/e m 

CTS 0.049 

0.8  -0.253 

1.0  -0.^54^ 

1.2  -0.589, 

1.3  -0-^25^ 
1.333  -0.634« 
2.0  -0.679 
2.74  -0.732^ 
4.0  -0.797« 
8.0              -0.813^ 

20.0  -0.768^ 


Prom  Reference  46. 


CHAPTBH  VII 
NUMERICAL  SOLUTIONS 

The  discussion  of  the  previous  chapter  has  made 
plausible  the  notion  that  the  bridge  diagram  set  B  may 
be  approximated  as  being  proportional  to  the  parallel 
diagram  set  P  ,  as  written  in  Eq.  (6.5).  The  PY  and  CHNC 
equations  are  then  but  two  particular  cases  of  this  approxi- 
mation in  which  m  is  always  taken  to  be  -1  and  0» 
respectively.  We  have  seen  how  to  allow  more  flexibility 
in  the  approximation  by  permitting  the  conditions  of  the 
problem  to  determine  m  ,  rather  than  assigning  it  a  value 
once  and  for  all.  The  arguments  presented,  however,  merely 
made  the  development  reasonable.   The  resulting  integral 
equations  are  notable  for  their  lack  of  "transparency"  and 
the  ultimate  test  of  any  of  them  lies  in  the  comparisons  of 
computed  answers  with  some  standard  solution. 

In  this  chapter  we  report  the  results  of  computations 
with  Eq.  (6.13),  which  we  will  call  the  "m"  equation,  for 
the  hard  sphere,  LJ  ,  and  Coulomb  potentials.  These 
computed  results  are  compared  in  each  case  to  the  answers 
obtained  from  the  PT  and  CHNC  equations,  as  reported  in  the 
literature  cited,  and  to  a  solution  taken  to  be  the  correct 
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answer.   The  latter  was  usually  a  Monte  Carlo  solution,  with 
the  exception  of  the  hard  sphere  equation  of  state,  as  noted 
helow.   Comparisons  with  experiment  are  not  reliable  tests 
of  the  accuracy  of  the  computed  results  so  long  as  the 
experimental  potential  is  not  accurately  known,  as  is  the 
present  case.  The  computed  functions  g  and  H  are  given 
in  tabular  form  in  Appendix  D. 

Before  discussing  the  results,  we  will  briefly  con- 
sider some  of  the  numerical  techniques  used  in  the  solutions. 

Once  a  specific  problem  has  been  decided  upon,  the 
parameter  m  must  be  found.  This  is  done  most  easily  in 
the  following  manner.   Define  the  function 

^(^12>  -  U(.r^^)fir^2^d^T^  (7.1) 

where  f  is  the  Player  function.   The  Fourier  transform  of 
this  equation  reads 

F(k)  -  f2(k)   .  (7.2) 

In  terms  of  P  ,  Eq,  (6.12)  becomes 


oo 


^(r)f(r)r^dr 


P^(r)r^dr 


m  -  -^B — »  (7.3) 

i 


2 
which  defines  the  average  of  f  over  the  distribution  F 

The  numerical  procedure  is  then  as  follows.   The  function  f 

is  computed  for  the  given  potential  and  temperature,  and  its 

Fourier  transform  is  obtained  using  the  techniques  of 


7^ 

Appendix  A.  Th.e  transform  of  P  Is  then  computed  from 
Bq.  (7.2)  and  the  Fourier  inversion  is  carried  out.   Finally 
Bq.  (7.3)  is  evaluated  using  Weddle's  rule  [^73.   An 
interesting  check  on  the  accuracy  of  the  numerical  procedure 
is  affoorded  by  the  Gaussian  f  function,  which  permits  m 
to  be  computed  analytically.   The  analytical  and  numerical 
computations  of  m  both  gave  the  same  number,  m  »  -0.35355» 
to  five  significant  figures.   This  accuracy  is  not  in  general 
to  be  expected,  however,  unless  the  interval  is  very  small, 
as  the  Gaussian  is  an  especially  smooth  function. 

The  "m"  equation  was  solved  in  its  Fourier  transform 
representation  by  an  iterative  scheme.  We  define  first 
the  quantities 

H(r)  .  g(r)e^<^(^)-l  (7.4) 

Q(r)  -  g(r)eP^^^^-  g(r)  (7.5) 

and 

P'(r)  -  (l+m)P(r)   .  (7.6) 

The  functions  Q,  p',  and  G  may  all  be  written  in  terms 
of  H  as  follows: 

Q(r)  .  g(r)eP<^(r)[l  -  e-Wr)^  (^.^^ 

-  -  f(r)[l  +  H(r)]  (7.8) 

p'(r)  -  (1+m)  [H(r)  -  ln[l  +  H(r)]]         (7-9) 


and 


G(r)  -  H(r)  -  Q(r)   .  (7.10) 
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Then  Eq.  (6.13)  may  "be  written  In  terms  of  H,  Q,  and  P  , 
giving 

H(r)  -  P'(r)  +  n  /  [P'(rO  -  Q(r')][H(|?  -  r'\   ) 

-  q(|r  -  t'\  )]d5r'   .         (7.11) 

Using  the  convolution  theorem  we  may  write  the  Fourier 
transform  representation  of  this  equation  as 

H(k)  -  P'(k)  +  n[P'(k)  -  Q  (k)][H  (k)  -  Q  (k)]   (7.12) 
whence  we  get  finally  that 

^     ^         [Q  (k)  -  P'(k)] 

H(k)  -  Q(k) ^ ~ .  (7.13) 

1  +  n[Q  (k)  -  P'(k)] 

Equations  (7.13) »  (7.8),  and  (7.9)  can  be  solved  by 
iteration.   The  procedure  is  as  follows.   An  initial  guess, 
H   ,  is  made  for  H  .   This  is  used  to  compute  Q   and  P 
from  Eqs.  (7.8)  and  (7.9).   The  Fourier  transforms  of  Q 
and  P   are  then  found  and  used  in  Eq.  (7.13)  to  obtain 
a  new  H  ,  denoted  H-,  .   This  quantity  H,   is  then  an 
input  quantity  for  the  computation  of  Q-,   and  P ,  ,  which 
in  turn  are  used  to  calculate  a  new  E^   »  and  so  on,  for 
any  desired  number  of  cycles. 

The  problem  is  considered  solved  when  a  solution  is 
self-consistent.   That  is,  when  the  new  H  calculated  is 
identical  to  the  H  used  as  input.   A  measure  of  the  degree 
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of  self-consistency  Is  given  by  the  root-mean-square  (rms) 
difference  between  successive  H's  ,  where  the  rms  is 
defined  as  follows: 


rms 


^  2  T  1/2 


(7.1^) 


Here  N  is  the  number  of  points  in  the  nximerical  solution, 
A  is  the  interval  between  points,  and  H,   is  the  value  of 
H  after  the  kti>  iteration.  In  practice,  a  solution  was 
considered  converged,  i.e.,  self-consistent,  if  the  rms 
between  successive  iterations  was  of  the  order  of  10  ^  or 
less.   An  additional  criterion  was  that  the  input  and  out- 
put H's  should  agree  to  at  least  two  significant  figures 
at  all  points. 

For  small  densities  the  rate  of  convergence  is  quite 
fast,  a  self-consistent  solution  usually  being  achieved  in 
about  10  iterations.   For  larger  densities,  the  slow  rate 
of  convergence  may  be  improved  by  the  following  device,  due 
to  Broyles  [A-8].  When  H  is  sufficiently  close  to  self- 
consistency,  it  is  found  to  approach  its  final  value 
exponentially.   Then  from  any  three  successive  iterations, 
k-1,  k,  and  k+1 ,   we  can  compute  the  value  that  H  would 
have  after  an  infinite  number  of  iterations,  viz., 


77 


H^l^CdA)  -  Hj^^^(jA)H^^^(dA) 


This  extrapolated  H  is  then  used  as  the  new  H  for  the 
input  of  the  next  iteration.   Some  care  must  be  taken  in 
employing  this  trick,  however.   If  it  is  used  too  soon,  i.e., 
before  the  solution  has  settled  down  to  an  approximately 
exponential  rate  of  approach,  it  will  cause  the  succeeding 
iterations  to  "blow  up." 

An  additional  feature  that  is  sometimes  useful  in 
speeding  convergence  is  the  so-called  mixing  process,  where 
we  write,  after  the  ktt  iteration, 

Hjj.''®*(jA)  -  a  Hj^(dA)  +  (l-a)Hj^^3^(dA)   .  (7.16) 

Normally  the  constant  a  is  set  equal  to  one.   VHien  a 
solution  tends  to  oscillate,  however,  a  smaller  a  will 
tend  to  dampen  the  oscillations  and  improve  the  convergence 
rate. 

The  Lennard-Jones  Potential 

The  LJ  potential  is  a  reasonably  realistic  potential 
which  still  retains  the  advantages  of  an  analytic  formula- 
tion.  It  is  given  by  the  expression 

Hr)    -  ^e[(a/r)^2  -  (a/r)^]  (7.17) 

where  the  constants  £   ,  having  the  dimensions  of  energy,  and 
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a  f  with  the  dimensions  of  length,  are  to  be  determined 
from  experiment.  Expressions  involving  the  LJ  potential 
may  be  put  in  dimensionless  form  by  taking  C     as  the  unit 
of  energy  and  a  as  the  unit  of  length.  Then  we  write 

30(r)  -  '<-(l/r^-l)/T*r^  (7.18) 

where  the  (dimensionless)  reduced  temperature  T*  is 
defined  by 

T*  -  kT/e  .  (7.19) 

The  corresponding  dimensionless  density  is  n*  •  na^  . 

Wood  and  Parker  [27]  have  used  the  Monte  Carlo  (MC) 
method  to  calculate  the  radial  distribution  function,  the 
pressure,  and  the  energy  for  a  fluid  interacting  with  a 
LJ  potential  at  several  densities  of  the  T*  ■  2.7^  iso- 
therm. Broyles  et  al.  [23 ♦'♦•9]  have  computed  the  solutions 
to  the  PT  and  GHNO  equations  at  these  same  densities  and 
temperature.  We  have  thus  chosen  these  cases  also  for 
purposes  of  comparison.   The  computed  value  of  m  was 
-0.732  for  this  case. 

The  results  of  the  "m"  equation  and  the  Monte  Carlo 
solution  are  shown  for  T*  «■  2.7^  and  n*  »  2/5,  5/6,  and 
1.0  in  Figures  7.1-7.5,  respectively.  The  agreement  is  in 
general  quite  good,  being  poorest  at  the  first  peak.   The 
PT  equation  also  yields  good  solutions  for  these  cases, 
while  the  OHNC  equation  is  in  poorer  agreement. 
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All  foxLT  plotted  curves  would  lie  confusingly  close 
and  they  have  been  compared  in  tabular  form  in  Tables  7»1- 
7.3.  These  tables  give  a  measure  of  the  deviation  from  the 
MO  solution  of  the  g's  computed  from  the  integral  equations 
in  the  form  (g^^^/g)-!  »  as  in  Reference  ^9.   All  three 
computed  g's  are  in  essential  agreement  at  the  reduced 
density  of  2/5  and  the  differences  in  Table  7.1  are  not  very 
meaningful.  More  substantial  deviations  appear  in  Table  7.2 
for  the  5/6  density  case,  as  well  as  for  the  1.0  density  in 
Table  7.5.   It  is  seen  that  the  present  equation  gives  the 
closest  agreement  in  general,  followed  in  turn  by  the  PT 
and  then  the  CHNC  equations. 

Once  the  radial  distribution  function  is  known  the 

pressure  p  and  energy  E  may  be  easily  found  from  the 

expressions  [1] 

00 
p/nkT  -  1  -  (27ni/5)  T 30'(r)g(r)r^dr       (7.20) 

and 

B/KkT  -  I  -  E'/NkT  (7.21) 

00 
-  2Tm  r  30(r)g(r)r^dr  .        (7.22) 

Values  of  pressure  and  temperature  for  the  three  oases 
mentioned  above  are  shown  in  Table  7.^.   The  MC  values  were 
chosen  from  the  various  computations  reported  by  Wood  and 
Peirker  [27]  on  the  basis  of  the  run  having  the  largest 
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TABLE  7.1 

COMPARISON  OP  PY,  CHNC,  AND  "m"  INTEGRAL  EQUATION  SOLUTIONS 
TO  MONTE  CARLO  FOR  THE  LJ  POTENTIAL:   T*  »  2.7^,  n*  =  2/5^ 


r 

(gpic/6m)-l 

(gj^Q/gpY)-^ 

^^MC^^CHNC^*^ 

0.90 

-0.027 

-0.016 

-0.07^ 

0.95 

-0.020 

-0.012 

-0.0^9 

1.00 

-0.066 

-0.060 

-0.079 

1.05 

-0.0^5 

-O.OA-1 

-0.0^8 

1.10 

-0.027 

-0.029 

-0.026 

1.15 

-0.011 

-0.013 

-0.006 

1.20 

0.003 

0.002 

-0.010 

1.3 

0.003 

0.001 

0.008 

1.^ 

0.006 

0.005 

0.007 

1.5 

0.001 

0.001 

-0.003 

1.6 

-0.009 

-0.006 

-0.013 

1.7 

-0.020 

-0.016 

-0.02^ 

1.8 

-0.012 

-0.013 

-0.01^ 

1.9 

-0.019 

-0.01^ 

-0.019 

2.0 

-0.012 

-0.012 

-0.012 

2.1 

-0.003 

-0.003 

-0.002 

2.2 

0.002 

0.000 

0.003 

^Tlie  PY  data  are  taken  from  Broyles,  Reference  49. 
The  CHNC  data  is  also  due  to  Broyles  (private  communication) 
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1.0 


Fig.  7. 2f— Radial  distribution  function  for  the  LJ  potential: 
T*  -  2.7^,  n*  =  5/6. 
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TABLE  7.2 

OCMPASISON  OF  PT»  CHNC,  AND  -m"  INTEGRAL  EQUATION  SOLUTIONS 
TO  MONTE  OAHLO  FOR  THE  LJ  POTENTIAL:   T*  -  2.?^,  n*  -  5/6* 


r 

^Wem)-^ 

(gj^Q/Epy)-! 

^^MC^^GHNC^"-^ 

0.90 

0.71^ 

0.772 

-0.117 

0.95 

0.055 

0.078 

-0.141 

1.00 

-0.089 

-0.083 

-0.138 

1.05 

-0.057 

-0.047 

-0.001 

1.10 

-0.021 

-0.041 

0.065 

1.15 

0.004 

-0.020 

0.102 

1.20 

0.007 

-0.016 

0.089 

1.3 

0.020 

0.020 

0.040 

1.^ 

0.007 

0.060 

0.007 

1.5 

0.011 

0.036 

-0.031 

1.6 

-0.015 

0.003 

-0.045 

1.7 

0.009 

0.014 

-0.010 

1.8 

0.010 

0.006 

-0.001 

1.9 

-0.004 

-0.007 

-0.009 

2.0 

-0.014 

-0.021 

-0.003 

2.1 

-0.002 

-0.010 

0.023 

2.2 

0.000 

0.004 

0.018 

^Ttie  PY  data  are  taken  from  Broyles,  Reference  49. 
The  CHNC  data  is  also  due  to  Broyles  (private  communication) 
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Fig.  7. 5. -Radial  distribution  function  for  the  LJ 
potential:   ?•  =  2.?^,  n*  =  1.0. 
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TABLE  7.5 

COnPARISON  OP  PY,  CHNC,  AND  "m"  INTEGRAL  EQUATION  SOLUTIONS 
TO  MONTE  CAHLO  POH  THE  LJ  POTENTIAL:   T*  =  2.?^,  n*  -  1.0^ 

r         (Sfic/8m>-^     ^WSPY^-^     (Smc/^Pj'^-^ 


0.90 

0.295 

0.555 

-0.504 

0.95 

-0.00^ 

0.024 

-0.214 

1.00 

-0.0^2 

-0.042 

-0.051 

1.05 

-0.068 

-0.085 

0.080 

1.10 

-0.050 

-0.081 

0.145 

1.15 

0.001 

-0.035 

0.169 

1.20 

0.047 

0.021 

0.157 

1.3 

0.102 

0.145 

0.051 

1.^ 

0.055 

0.116 

-0.082 

1.5 

-0.009 

0.025 

-0.097 

1.6 

-0.057 

-0.042 

-0.080 

1.7 

-0.012 

-0.050 

-0.025 

1.8 

-0.015 

-0.051 

-0.017    . 

1.9 

-0.024 

-0.055 

-0.012 

2.0 

-0.007 

-0.014 

0.051 

2.1 

0.000 

-0.002 

0.055 

2.2 

0.045 

0.065 

0.044 

*The  PT  data  are  taken  from  Broyles,  Reference  49. 
The  CHNC  data  is  also  due  to  Broylea  (private  communication) 
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TABLE  7.^ 
THERMODYITAMIG  QUANTITIES  FOR  THE     LJ     POTENTIAL 


n* 

2/5 

5/6 

1.0 

(p/nlcT)j^Qa 

1.2-1.5 

^.08 

6.97 

(p/nkT)^ 

1.2^ 

^.09 

6.99 

(p/nkT)pYb 

1.2^ 

4.01 

6.8 

(p/nkT)(.gjjcb 

1.28 

5.11 

9.1 

(E;'NkT)„Qa 

-0.86 

-1.57 

-1.60 

il&'/WLT)^ 

-0.865 

-1.58 

-1.615 

(K>NkT)pYt 

-0.865 

-1.61 

-1.67 

(E;^NkT)(,gjjcb 

-0.859 

-l.'l-O 

-1.19 

^rom  Reference 

27. 
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number  of  moves  per  particle,  which  might  be  expected  to 
give  the  best  average.  Here  again  the  present  equation 
gives  the  closest  answers,  though  the  differences  from  PY 
are  small. 

The  results  for  the  "m"  equation  in  Table  7.^  include 
a  small  correction  for  the  truncation  of  the  integrals  in 
Eqs.  (7.20)  and  (7.22).   This  truncation  is  unavoidable  in 
machine  calculations  and  was  corrected  for  in  the  following 
way.   If  the  maximum  value  of  r  in  the  numerical  solution 
is  R  ,  we  write,  e.g., 

R 

EVnkT  -  2iin  /       30(r)g(r)r^dr 
■Jo 

CO 

30(r)g(r)r2dr  (7.23) 

H 


+  2Tm  J 


-   (27iikT)^^^p  *  Ag     ,  (7.2^) 

Here  the  first  term  on  the  right  is  the  value  obtained  on 
the  computer.  The  correction, 

Ag  -  2Tm  I   00(r)g(r)r''dr  ,  (7.25) 

may  be  approximated  by  setting  g  equal  to  one  everywhere 
In  the  region.   If  R  is  sufficiently  large  this  approxi- 
mation will  be  quite  good.   For  the  LJ  potential,  then, 
we  easily  obtain 


^B 

as   — 

8Tm' 
5T»R^ 

(1 

-  1/3H^) 

Sto  - 

87m* 

» 

88 
(7.26) 

(7.27) 


in  reduced  units,   A  similaLT  correction  was  applied  to  the 
machine  solutions  for  the  pressure.  For  the  n*  «  1  case, 
for  example,  we  took  R  ■  ^.5  (in  units  of  a).  This  gave 
a  computed  sVnkT  of  -1.582  and  a  correction  A™  of 
-0,033,  or  about  2  per  cent  of  the  total. 

The  Coulomb  Potential 

We  consider  here,  following  Carley  C50,51]»  a  system 
of  classical  point  charges  embedded  in  a  rigid  uniform 
background  of  charge  of  the  opposite  sign,  such  that  the 
net  charge  is  zero.  The  pair  potential  is  Just  the  Coulomb 
potential 

0(r)  -  e^/r  (7.28) 

where  e  is  the  charge  on  an  electron.   For  convenience  we 
will  use  the  same  reduced  units  employed  in  the  above 
references,  viz.,  the  unit  of  length  is  taken  to  be  the 
ion  sphere  radius  a  given  by 

a  -  [3/(^rai)3^/^  (7.29) 

and  the  dimensionless  reduced  "temperature"  is 


©  -  kTa/e^   ,  (7.30) 


so  that  we  may  write 
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P0(r)  -  l/(©r)  (7.51) 

where  r  is  in  units  of  a  .   Note  that  ©  is  a  function, 
of  the  density  as  well  as  of  the  temperature. 

Equation  (6.15)  has  been  solved  for  the  Coulomb 
potential  for  the  two  cases  ©  -  O,'*-  and  ©  ■  1.0.   The 
computed  values  of  m  for  these  cases  were  -0.566  and 
-O.25I,  respectively.   It  was  pointed  out  in  the  last 
chapter  that  due  to  the  divergent  nature  of  integrals  con- 
taining Coulomb  f  functions,  the  corresponding  Debye- 
Htlckel,  or  exponentially  shielded  Coulomb,  potentials  were 
used  to  compute  m  . 

The  computed  results  are  given  in  Figures  7.^  and 
7.5»  where  comparison  is  made  to  solutions  of  the  PY  said 
OHNC  equations  obtained  by  Carley  [50,51]  and  to  Monte 
Carlo  calculations  by  the  same  author  [26],  A  clear  improve- 
ment over  the  PT  and  CHNC  results  is  apparent  in  both  cases 
by  use  of  the  self-adjusting  m-approximation  of  the  last 
chapter.   The  results  of  the  ©  -  1.0  case  axe  especially 
good,  the  computed  g  diverging  from  the  Monte  Carlo 
average  only  in  the  region  where  g  approaches  one. 

The  Hard  Sphere  Potential 

As  a  physical  model,  a  fluid  of  hard  spheres  is 
adequate  only  for  describing  systems  at  very  high  tempera- 
tures.  Nonetheless,  the  simple  analytical  form  of  the  hard 
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sphere  potential  lias  invited  ei^austive  studies  of  this 

model,  in  that  otherwise  intractable  expressions  take  on 

manageable  forms. 

For  hard  spheres  of  diameter  a  ,  the  potential  is 

given  by 

I  oo  »   r  <  a 
30(r)  -  {  (7.32) 

I  0   ,   r  >  a 

and  is  temperature  independent.  We  will  use  the  hard  sphere 

diameter  a  as  the  natural  unit  of  length,  so  that  the 

reduced  density  is  n*  -  na^  .  Prom  Bq.  (7.32)  it  is  easy 

to  see  that  the  Mayer  f  function  is  here  a  step  function, 

viz., 

-1   ,   r  <  1 
f(r)  -  \  (7.33) 

0   ,   r  >1 

where  r  is  in  units  of  a  .  The  derivative  of  f  is 
thus  the  Dirac  delta  function, 

f'(r)  -6(r-l)  (7.5^) 

and  80,  from  Eq.  (C.7)  we  get  for  the  (virial)  equation  of 
state 

j||  -  1  ^  ^  n*g(l)   .  (7.35) 

We  have  solved  Eq.  (6.15)  for  the  hard  sphere  model 
for  nine  densities,  from  n*  -  0.1  to  0.9  ,  and  have 
obtained  the  pressure  in  each  case  using  Eq.  (7.55).  The 
computed  g's  and  H's  are  given  in  Tables  D.5  to  D.7.   The 
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computed  equation  of  state  is  shown  in  Pig.  7.6,  along  with 
the  PT  and  CHNC  equations  of  state  obtained  "by  Klein  [24, 
25],  also  using  Bq.  (7.55).  The  reference  isotherm  was  put 
together  by  Klein  by  taking  the  equation  of  state  obtained 
from  the  five  known  virial  coefficients  for  n*  less  than 
0.5  and  the  molecular  dynamics  calculations  for  n*  greater 
than  0.7.  For  0.5  ^  n*  ^  0.7,  the  reference  isotherm  is  a 
graphical  extrapolation  of  the  former  onto  the  latter. 

Numerical  solutions  of  the  integral  equations  for  g 
for  the  hard  sphere  potential  are  extremely  sensitive  to  the 
size  of  the  interval,  as  has  been  pointed  out  by  Sahlin 
C56].  It  is  found  in  practice  that  g  at  the  point  1 
will  decrease  as  the  interval  size  is  decreased.   Thus  the 
solution  shown  in  Pig.  7.6,  computed  for  an  interval  size 
of  0.025  (in  reduced  units),  probably  lies  slightly  above 
what  would  be  the  correct  solution  for  an  infinitesimal 
interval  size.  Nonetheless,  it  is  possible  to  conclude,  on 
the  basis  of  the  available  numbers,  that  the  •'m"  equation 
yields  improved  solutions  over  the  PY  and  CHNC  equations 
for  the  hajxi  sphere  potential. 


9^ 


0  0.^  0.8      ^,         1.2 

Fig.  7. 6. -The  equation  of  state  of  a  fluid  of  hard  spheres. 


APPENDICES 


APPENDIX  A 

FOURIER  TRANSFORMS 

There  is  a  certain  amount  of  ambiguity  in  the 
definition  of  Fourier  transforms  in  that  factors  of  2n 
may  be  distributed  in  several  ways.   This  is  harmless  so 
long  as  a  given  definition  is  applied  consistently.  We 
have  used  throughout  this  dissertation  the  following 
definitions: 


F(r) 


P(k) 


(2ti) 


^  rP(k)e^^*^d^k  (A.l) 


jF(r)e"^^*Vr   .  (A. 2) 


If  F  is  a  function  of  the  magnitude  of  r  alone, 

then  F  is  a  function  only  of  the  magnitude  of  k  ,  and  the 

integrations  over  angles  above  may  be  carried  out,  giving 

oo 

p(r)  .  — i-  /  F(k)  k  sin  krdk  (A. 3) 

P(1j:)  »  ^   /  F(r)  r  sin  krdr   .  (A.A-) 

i 

The  Fourier  transform  of  the  Coulomb  l/r  potential 
may  be  found  most  easily  by  means  of  a  well-known  trick. 
Write 
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l/r  -     11m     e'^/r     .  (A.  5) 

a  -»  o 

Then  the  transform  of  l/r  may  be  written 

lim     ^    /  e'^  sin  krdr  -     lim     ^  -^^  (A. 6) 

a-ooj  a-o-o         a+k 

0 

-  ^hA^      .  (A. 7) 

The  Fourier  transform  of  30(r)   for  the  Coulomb  potential 
is  thus 

30(k)  -  ^n3e^/k^  ,  (A. 8) 

as  stated  in  Chapter  II. 

For  numerical  work  the  infinite  sum  over  infinitesimal 
intervals,  represented  by  the  symbol  /  dx  ,  must  be  replaced 
by  a  finite  sum  over  finite  intervals.   So  long  as  only  the 
value  of  an  integral  is  desired,  the  machine  calculation 
may  be  performed  by  the  best  convenient  rule,  such  as 
Weddle's  rule  [^73.   For  Fourier  integrals,  however,  more 
care  must  be  taken.   The  basic  property  of  the  Fourier 
transform  that  makes  it  eminently  useful  is  the  fact  that 
it  is  an  expansion  in  orthogonal  functions,  making  possible 
the  reciprocal  relation  of  a  function  to  its  transform 
indicated  by  Eqs.  (A. 5)  and  (A.^).   If,  in  a  given  numerical 
computation,  this  reciprocal  property  is  to  be  used,  the 
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discrete  representation,  of  the  expansions  must  maintain 

orthogonality.  We  will  now  consider  how  this  is  done.   The 

following  technique  is  due  to  Professor  A.  A.  Broyles. 

Choose  a  range  (0,r   )  for  the  r  variable  and 
^  '  max 

(0,k   )  for  the  k  variable.  Placing  N  equally-spaced 
points  in  the  range  of  x  divides  it  into  N  intervals  of 
equal  size  Ar  »  r   /N  .  A  continuous  function  F(r)   is 
then  represented  by  the  discrete  function  F(nAr)  ,  n  =  0,N 
Similarly,  the  range  of  k  is  divided  into  M  equal 
intervals  dJs.   o  ^max^^  *   Equations  (A. 5)  and  (A.^)  then 
become,  respectively,  by  the  trapezoidal  rule  [A-7] , 

2    " 
P(nAr)  -  ^^^ —  y      F(mAk)m  sin  (mnArAk)      (A. 9) 

2       ^ 
F'(niAk)  -  ^^^^^    ^   F(nAr)n  sin  (mnArAk)  .   (A. 10) 

n»l 

(The  m  »  0  and  n  »  0  terms  are  not  written  in  the  sums 
since  they  are  identically  zero.)  Thus  we  have 

M 


p(n^)  .  2^Ar  ^   F(n'Ar)n'  Y       sin  (mn' ArAk 
n'  »1  m^ 


X  sin  (mnArAk)  .  (A. 11) 

Consider  now  the  second  sum  in  Eq.  (A. 11), 


J  •  >    sin  (mn'ArAk)  sin  (uinArAk) 
m-1 
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(A. 12) 


M 


-  i   V  cos[(n' -n)mArAk]  -  cos[(n' +n)mArAk] 


m=l 
With  the  identity  [52] 

M 

2   ^   cos  m©  »  -  1  + 
m«l 

we  may  write 


sin  (M  +  ^)© 

T f 

sin  I  © 


(A. 15) 


(A. 14) 


J  - 


sin  I  (2M+l)(n'-n)ArAk 
sin  -H   (n'-n)ArAk 

sin  I  (2M+l)(n'+n)ArAk 
sin  ^  (n' +n)ArAk 


(A.15) 


We  now  want  to  impose  orthogonality.  Write 

J  -  c5„„'  (A. 16) 


nn 


where  U aa      ^^  "^^^   Kronecker  delta, 


5u- 


(A. 17) 


1   ,   i  "  J 

0   ,   i  /  d 

and  C  is  a  constant  to  be  determined.   Then  for  n  /  n% 
we  require 
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sin  i  (2M+l)ArAk(n'-ii)   sin  i  (2M+l)ArAk(n' +a) 

2   .  ^~^r ^ .(A. 18) 

sin  4  ArAk(n'-n)         sin  ^  ArAk(n' +n) 

This  is  satisfied  for  arbitrary  integers  n  and  n'  , 
n  ^  n'  ,  if 

(2M+l)ArAk  -  2n  (A. 19) 

since  both  sides  then  vanish  identically.  For  n'  =  n,  the 
second  term  in  Eq.  (A. 15)  vanishes  also,  using  the  result 
Bq.  (A. 19).   The  first  term  gives 


n  -.►n'  sin  ^ 


li^  si^  ;  i^-^'A    ,   lim    n  cos  n   (^-O        (A.20) 
1'  sin  ^^^         n^n'  ^   cos  .  i|^ 

-  2M  +  1  .  (A.21) 


Thus, 

and  Eq.  (A. 11)  becomes 

nn^)   .  ^^  ^    ^     n  Fin  ^)6^.  (A.25) 

n'  -1 

-  P(nAr)  .  (A. 2^) 

If  we  had  started  from  P  instead  of  F  ,  we  would 
have  found  that  we  needed  the  condition 

(2N+l)ArAk  -  2n  '  (A.25) 
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for  orthogonality.   Since  the  intervals  are  the  same,  this 
implies  ve  must  have  N  »  M  . 

Summing  up,  then,  we  find  that  the  reciprocal  set  of 
Fourier  transforms  for  continuous  functions,  Bqs.  (A. 3)  ancL 
(A.^),  is  adequately  represented  in  discrete  form  by  the 
following  set: 

2    ^ 
P(n)  -  — A —  y    m  P  (m)  sin  (nmA)  (A. 26) 

m>l 

N 


^  ^"'^  ■  ^  ^  n  P  (n)  sin  (nmA)  (A.27) 

mal 


where 


and 


A  -  27i/(2N+l)  (A. 28) 


D  =  (Ar)^   .  (A. 29) 

Note  that  once  the  range  and  interval  size  of  the  r 
variable  are  chosen,  the  orthogonality  condition  fixes  the 
range  and  interval  size  of  the  k  variable.   Furthermore, 
the  development  above  can  only  be  carried  out  for  the 
trapezoidal  rule  of  numerical  integration,  which  gives  equal 
weights  to  all  the  points  [^7]. 


APPENDIX  B 

PARAMETERS  IN  THE  LEBOWITZ  SOLUTION  FOR  A  MIXTURE  OP 

HARD  SPHERES 

Lebowitz's  exact  solution  of  the  PY  equation  for  a 
mixture  of  hard  spheres  [^0]  gives  the  three  direct  cor- 
relation functions  of  the  system,  Eq.s.  (5.7)  and  (5.8), 
as  a  function  of  r  and  seven  parameters,  called  here 
a,b,a,b,b,d,  and  \  .   These  parameter  are  given 
below. 

Let  n   and  n   be  the  number  densities  of  the 

+       - 

positive  and  negative  hard  sphere  ions,  respectively,  and 

R   and  R   the  corresponding  hard  sphere  diameters.  Ve 

will  assume  that  R  ^  R   .   We  then  define  the  following 
-  —  + 


terms: 


and 


w  -  Ttn  /6  (B.l) 

+     + 


\  -  (R_  -  R^)/2  (B.2) 


R^_-  (R_  +  R_^)/2  (B.3) 


z-wR^+wR^   .  (B.^) 

+  +     -  - 
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Lebowitz  found  that  the  pressure  p  of  the  mixture  la 
given,  in  the  PY  approximation,  by 


0P 


[  ^%  *  ^- 


)(1  +  z  +  z'^) 


i-  w  w  (R  -  R  )^[R  +R  -(-RR(wR^ 


'1/ 


^.  w_R  2)]  I  /  (1  _  2)3 


(B.5) 


The  hard  core  radial  distribution  functions,  labelled 
here  with  a  "short-range"  superscript,  were  found  to  have 
the  following  values  at  the  hard  core  diameters: 

6^^(R^)  -  [1  +  ^  z  ^-  I  w_R_2(R^  .  R_)]  /(I  -  2)2    (B.6) 

gffCRj  -  [1  +  I  z  +  I  w^R^2(g^  .  R^)]  /  (1  _  2)2    (B,7) 


s!!(V)  -  fR-s+H.(R+)  +  V— ^^-^^/^V    • 


(B.8) 


With  these  defined  quantities,  we  complete  the  set 
of  required  parameters  by  listing  the  following  quantities 


1  +  z  +  z^  +  g  \^(^+  *   n_)(l  +  2z) 


-  3w  (R   -  R  )2[R   +  R   +  R  R  (2w  R  2 


■••  w  R  ^)] 


/  (1  -  z)5  ^.  ^  R^^0p/(1  -  z) 


(B.9) 


10* 

b^  -  -  6Cw  R  2g  2^^j  ^  ^  ^     2   2^   j^  ^3^^Q^ 


b  -  -  6[w  R  g_(R  )  +  w  R  g  (R  )]  (B.ll) 

d  -  |[w^a_^  +  w_a_]  .  (B.12) 

Tbe  constants  a   and  b   axe  obtained  from  a   and  b 
by  interchanging  w  ,R   with  w_,R_  . 


APPENDIX  C 
THE  VIRIAL  AND  COMPRESSIBILITY  EQUATIONS  OP  STATE 

The  equation  of  state  of  a  fluid  composed  of  molecules 
which  interact  through  spherically  symmetric  pair  forces  may 
be  obtained  in  two  ways,  one  proceeding  from  the  virial 
theorem,  the  other  based  on  the  fluctuations  in  density. 
Derivations  of  both  expressions  may  be  found  in,  e.g., 
Hirschfelder,  Curtiss,  and  Bird  [53].   They  are  (a)  from 
the  virial  theorem, 

^  .  ,  .  2^  f  3  d|lrl  g(^)^3^^  (,^,) 

and  (b),  from  the  compressibility, 

kT(dp/^n)^^  -  1  +  ^  Tin  /  [s(r)-l]r^dr  .         (0.2) 

Equation  (C.2)  may  be  put  in  a  more  convenient  form 
by  using  the  fact  that  the  Fourier  transforms  of  G  and  T 
are  related  by  the  expression  [see  Chap.  I] 

1  +  nG(k)  -  1/[1  -  nT(k)]   ,  (C.3) 

so  that 

kTCdp/^n)"-'-  -  1  +  nG(0)  (C.^) 

-  1/[1  -  nT(0)]  (C.5) 
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and 


^  ^  -  1  -  ^  Ttn  /  T(r)r^(ir  .  (0.6) 

J   A 


0 

Equation  (C.l)  may  also  be  conveniently  rewritten  as  follows, 

-^  =  1.2^/  g(r)eP^^^)f'(r)r5dr  (C.7) 

where  f '  is  the  derivative  of  the  Mayer  f  function. 

It  can,  of  course,  be  shown  [5^]  that  if  the  exact 
diagram  expansions  of  T  and  ge''^   are  inserted  in  the 
equations  (0.6)  and  (0.7) »  both  give  the  correct  virial 
expansion.  When  approximations  to  the  bridge  set  are  intro- 
duced, however,  the  two  equations  give  inconsistent  answers. 
This  situation  could  be  exploited  to  determine  the  parameter 
m  of  Chapter  V  by  imposing  consistency  between  Eqs.  (0.6) 
and  (0.7)  even  with  the  approximated  B  set.   That  is,  we 
write 

/  T(r)r^dr  '  -\\   g(r)e^'^^^^f '(r)r^dr 

-  ^  n  ^  rg(r)e^^^^)f'(r)r5dr  .   (0.8) 

Using  the  diagram  expansions  of  Eqs.  (I.30)  and  (1.29), 
we  find  that  this  leads  to  the  condition 


L 
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oo 

[p(r)  +  B(r)  +  f(r)[P(r)+B(r)]  ^  |  f '(r)CP(r)+B(r)] 
^[P(r)H.B(r)]l  Ar  -  -  /  [f(r)S(r) 


>|  f'(r)  n  ^ 


0 

+  I  f  (r)S(r)  +  I  f '(r)  n  ^isilr^dr  .      (C.9) 

If  B  is  now  approximated  by  mP  ,  we  get  for  m  the 

equation 

oo 

(1+m)  /  [P(r)  +  f(r)P(r)  +  |  f  (r)P(r) 

oo 
*  I  ^'(^^  ^T^^""^^  "  "  /  t^^^^^S^^)  *  f  f'(r)S(r) 

^  I  f'(r)  n  ^^Ir^dr  .  (CIO) 

The  sets  P  and  3  may  then  be  related  to  g  obtained 
from,  say,  the  PY  equation  and  the  integrals  evaluated  to 
give  m  . 


APPENDIX  D 

THE  CMIPUTBD  FUNCTIONS  g  AND  H  FOR  THB 
LJ,  COULOMB,  AND  HARD  SPHERE  POTENTIALS 

In  this  appendix  we  present  the  results  of  compu- 
tations for  the  radial  distribution  function  g  and  the 
function  H  defined  in  Bq.  C?.'*-)  from  the  integral 
equation  presented  in  Chapter  VI,  Eq,.  (6.13).  Table  D.l 
lists  the  three  LJ  cases  studied.   The  Coulomb  solutions 
are  given  in  Table  D.2,  and  the  hsurd  sphere  solutions 
take  up  the  remaining  tables,  D.3  to  D.7. 
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TABLE  D.l 
THE  COMPUTED  FUNCTIONS  g  AND  H  FOR  THE  LJ  POTENTIAL 


"  ""  ■ 

1^ 

2.7^ 

'f  *'   - 

2,7^ 

T^ 

2.74 

n*   - 

2/5 

n*   - 

5/6 

n*   - 

1.0 

r/a 

g 

H 

g 

H 

g 

H 

0.05 

0.000 

3.351 

0.000 

20.129 

0.000 

40.525 

0.10 

0.000 

3.221 

0.000 

19.285 

0.000 

58.558 

0.15 

0.000 

3.035 

0.000 

18.105 

0.000 

56.080 

0.20 

0.000 

2.819 

0.000 

16.769 

0.000 

55.5^1 

0.25 

0.000 

2.593 

0.000 

15.392 

0.000 

50.548 

0.30 

0.000 

2.568 

0.000 

14.052 

0.000 

27.804 

0.35 

0.000 

2.150 

0.000 

12.715 

0.000 

25.147 

0.^ 

0.000 

1.941 

0.000 

11.451 

0.000 

22.587 

0.^5 

0.000 

1.743 

0.000 

10.245 

0.000 

20.129 

0.50 

0.000 

1.557 

0.000 

9.092 

0.000 

17.776 

0.55 

0.000 

1.382 

0.000 

7.999 

0.000 

15.552 

0.60 

0.000 

1.128 

0.000 

6.964 

0.000 

15.405 

0.65 

0.000 

1.065 

0.000 

5.990 

0.000 

11.404 

0.70 

0.000 

0.922 

0.000 

5.079 

0.000 

9.557 

0.75 

0.000 

0.789 

0.000 

^.255 

0.000 

7.815 

0.80 

0.000 

0.666 

0.000 

5.457 

0.000 

6.243 

0.85 

0.003 

0.553 

0.006 

2.755 

0.010 

4.856 

0.90 

0.129 

0.450 

0.277 

2.126 

0.408 

5.601 

0.95 

0.663 

0.356 

1.260 

1.578 

1.755 

2.546 

1.00 

1.272 

0.272 

2.111 

1.111 

2.675 

1.675 

1.05 

1.578 

0.197 

2.275 

0.725 

2.615 

0.984 

1.10 

1.619 

0.131 

2.021 

0.412 

2.095 

0.462 

1.15 

1.536 

0.077 

1.672 

0.169 

1.557 

0.088 

1.20 

1.418 

0.024 

1.565 

-0.014 

1.158 

-0.164 

1.25 

1.303 

-0.018 

1.155 

-0.146 

0.899 

-0.522 

1.30 

1.205 

-0.052 

0.970 

-0.257 

0.747 

-0.412 

1.35 

1.125 

-0.080 

0.862 

-0.295 

0.667 

-0.454 

1,W 

1.063 

-0.101 

0.797 

-0.527 

0.656 

-0.462 

1.^5 

1.015 

-0.117 

0.762 

-0.558 

0.658 

-0.445 

1.50 

0.980 

-0.128 

0.751 

-0.552 

0.665 

-0.408 

1.60 

0.938 

-0.156 

0.780 

-0.282 

0.767 

-0.295 

1.70 

0.927 

-0.126 

0.856 

-0.192 

0.910 

-0.141 

1.80 

0.938 

-0.100 

0.966 

-0.074 

1.077 

0.053 

1.90 

0.968 

-0.061 

1.089 

0.056 

1.255 

0.196 

2.00 

1.002 

-0.021 

1.156 

0.150 

1.247 

0.219 

2.10 

1.025 

0.006 

1.122 

0.104 

1.100 

0.082 

2.20 

1.027 

0.014 

1.041 

0.027 

0.947 

-0.065 

2.30 

1.021 

0.011 

0.972 

-0.057 

0.875 

-0.155 

2.40 

1.012 

0.004 

0.939 

-0.068 

0.875 

-0.151 
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Table  D.l   (continued) 


T*   - 

2.74 

T*    . 

2.74 

-f*    . 

2.74 

n*    - 

2/5 

n*   - 

5/6 

n*    - 

1.0 

r/a 

S 

H 

S 

H 

g 

H 

2.50 

1.003 

-0.003 

0.939 

-0.067 

0.922 

-0.083 

2.60 

0.997 

-0.008 

0.960 

-0.0^5 

0.986 

-0.018 

2.70 

0.99^ 

-0.009 

0.989 

-0.01^ 

l.O^il- 

0.0^0 

2.80 

0.995 

-0.008 

1.016 

0.013 

1.076 

0.072 

2.90 

0.997 

-0.005 

1.051 

0.028 

1.071 

0.068 

3.00 

1.000 

-0.002 

1.029 

0.027 

1.035 

0.031 

3.50 

1.001 

0.000 

0.992 

-0.009 

0.993 

-0.008 
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TABLE  D.2 
THE  COMPUTED  FUNCTIONS  g  AND  H  FOR  THE  COULOMB  POTENTIAL 


©   - 

0.4 

6   » 

1.0 

r/a 

E 

H 

S 

H 

0.10 

0.000 

5.094 

0.000 

15.697 

0,20 

0.026 

2.820 

0.000 

14.527 

0.50 

0.127 

2.568 

0.005 

15.047 

0.^0 

0.274 

2.555 

0.025 

11.856 

0.50 

0.422 

2.121 

0.079 

10.755 

0.60 

0.555 

1.927 

0.167 

9.744 

0.70 

0.659 

1.752 

0.276 

8.821 

0.80 

0.745 

1.595 

0.595 

7.982 

0.90 

0.808 

1.454 

0.511 

7.221 

1.00 

0.856 

1.528 

0.618 

6.55^ 

1,10 

0.895 

1.217 

0.712 

5.913 

1.20 

0.920 

1.118 

0.791 

5.555 

1.30 

0.9^1 

1.050 

0.855 

4.850 

1.^0 

0.956 

0.955 

0.905 

^.597 

1.50 

0.967 

0.884 

0.945 

5.991 

1.60 

0.976 

0.825 

0.970 

5.628 

1.70 

0.982 

0.768 

0.989 

5.505 

1.80 

0.987 

0.720 

1.001 

5.014 

1.90 

0.990 

0.677 

1.008 

2.756 

2.00 

0.995 

0.658 

1.011 

2.527 

2.40 

0.999 

O.5I6 

1.004 

1.844 

2.80 

1.001 

0.450 

0.995 

1.451 

5.20 

1.000 

0.560 

1.003 

1.195 

5.60 

1.000 

0.287 

1.054 

1.071 

4.00 

1.000 

0.186 

1.047 

0.957 

4.50 

1.000 

0.009 

1.000 

O.6O5 
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TABLE  D.3 

THE  COMPUTED  FUNCTIONS  g  AND  H  FOR  THE  HARD  SPHERE 

POTENTIAL 


n* 

-  0.1 

n* 

-  0.2 

r/a 

S 

H 

6 

H 

0.1 

0.000 

0.525 

0.000 

1.379 

0.2 

0.000 

0.476 

0.000 

1.243 

0.3 

0.000 

0.431 

0.000 

1.112 

0.4 

0.000 

0.386 

0.000 

0.984 

0.5 

0.000 

0.343 

0.000 

0.862 

0i6 

0.000 

0.301 

0.000 

■    0.745 

0.7 

0.000 

0.261 

0.000 

0.634 

0.8 

0.000 

0.223 

0.000 

0.531 

0.9 

0.000 

0.188 

0.000 

0.436 

1.0 

1.155 

0.155 

1.349 

0.349 

1.1 

1.125 

0,125 

1.271 

0.271 

1.2 

1.098 

0.098 

1.202 

0.202 

1.5 

1.075 

0.073 

1.144 

0.144 

1.4 

1.052 

0.052 

1.094 

0.094 

1.5 

1.054 

0.034 

1.054 

0.054 

1.6 

1.020 

0.020 

1.023 

0.023 

1.7 

1.009 

0.009 

1.001 

0.001 

1.8 

1.001 

0.001 

0.987 

-0.013 

1.9 

0.996 

-0.004 

0.982 

-0.018 

2.0 

0.996 

-0.004 

0.985 

-0.015 

2.1 

0.997 

-0.003 

0.992 

-0.008 

2.2 

0.999 

-0.001 

0.997 

-0.003 
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TABLE  D.^ 

THE  COMPUTED  FUNCTIONS  g  AND  H  FOR  THE  HARD  SPHERE 

POTENTIAL 


n* 

-  0.5 

n* 

.   0.4 

r/a 

g 

H 

g 

H 

0.1 

6,666 

2.789 

0.000 

5.086 

0.2 

0.000 

2.494 

0.000 

4.512 

0.5 

0.000 

2.207 

0.000 

5.958 

0.4 

0.000 

1.952 

0.000 

5.427 

0.5 

0.000 

1.668 

0.000 

2.921 

0.6 

0.000 

1.418 

0.000 

2.445 

0.7 

0.000 

1.185 

0.000 

2.002 

0.8 

0.000 

0.969 

0.000 

1.597 

0.9 

0.000 

0.772 

0.000 

1.255 

1.0 

1.596 

0.596 

1.915 

0.915 

1.1 

1.445 

0.445 

1.645 

0.645 

1.2 

1.51^ 

0.514 

1.425 

0.425 

1.5 

1.206 

0.206 

1.255 

0.253 

1.4 

1.120 

0.120 

1.122 

0.122 

1.5 

1.055 

0.055 

1.029 

0.029 

1.6 

1.007 

0.007 

0.967 

-0.052 

1.7 

0.976 

-0.024 

0.954 

-0.066 

1.8 

0.961 

-0.059 

0.924 

-0.076 

1.9 

0.959 

-0.041 

0.955 

-0.065 

2.0 

0.972 

-0.028 

0.966 

-0.054 

2.1 

0.991 

-0.009 

1.005 

0.003 

2.2 

1.002 

0.002 

1.019 

0.019 

2.5 

1.007 

0.007 

1.022 

0.022 

2.4 

1.008 

0.008 

1.018 

0.018 

2.6 

1.004 

0.004 

1.004 

0.004 

2.8 

1.000 

0.000 

0.995 

-0.005 

11-* 


TABLE  D.5 

THE  COMPUTED  FDUCTIONS  g  AND  H  FOR  THE  HARD  SPHERE 

POTENTIAL 


n* 

-  0.5 

n* 

-  0.6 

r/a 

6 

H 

S 

H 

0.1 

0.000 

S.^21 

0.000 

16.561 

0.2 

0.000 

7.257 

0.000 

12.961 

0.3 

0.000 

6.296 

0.000 

11.195 

0.4 

0.000 

5.620 

0.000 

9.512 

0.5 

0.000 

4.552 

0.000 

7.923 

0.6 

0.000 

5.758 

0.000 

6.456 

0.7 

0.000 

5.024 

0.000 

5.066 

0.8 

0.000 

2.557 

0.000 

5.852 

0.9 

0.000 

1.767 

0.000 

2.755 

1.0 

2.259 

1.259 

2.852 

1.852 

1.1 

1.842 

0.842 

2.159 

1.159 

1.2 

1.515 

0.515 

1.612 

0.612 

1.3 

1.270 

0.270 

1.250 

0.250 

1.4 

1.096 

0.096 

1.019 

0.019 

1.5 

0.982 

-0.018 

0.887 

-0.115 

1.6 

0.915 

-0.085 

0.828 

-0.172 

1.7 

0.886 

-0.114 

0.825 

-0.177 

1.8 

0.890 

-0.110 

0.858 

-0.142 

1.9 

0.920 

-0.080 

0.924 

-0.076 

2.0 

0.974 

-0.026 

1.016 

0.016 

2.1 

1.050 

0.050 

1.095 

0.095 

2.2 

1.047 

0.047 

1.096 

0.096 

2.3 

1.042 

0.042 

1.061 

0.061 

2.4 

1.026 

0.026 

1.020 

0.020 

2.5 

1.009 

0.009 

0.987 

-0.015 

2.6 

0.996 

-0.004 

0.969 

-0.051 

2.8 

0.985 

-0.015 

0.972 

-0.028 

5.0 

0.994 

-0.006 

1.002 

0.002 

5.2 

1.004 

0.004 

1.017 

0.017 
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TABLE  D.6 

THE  COMPUTED  FUNCTIONS  g  AND  H  FOR  THE  HARD  SPHERE 

POTENTIAL 


n* 

-  0.7 

n* 

-  0.8 

r/a 

S 

H 

g 

H 

0.1 

C>M6 

2^.0^7 

6.000 

55.622 

0.2 

0.000 

22.022 

0.000 

28.156 

0.3 

0.000 

17.170 

0.000 

25.995 

OA 

0.000 

1^.^77 

0.000 

20.096 

0.5 

0.000 

11.9^0 

0.000 

16.444 

0.6 

0.000 

9.575 

0.000 

15.052 

0.7 

0.000 

7.^07 

0.000 

9.961 

0.8 

0.000 

5.^70 

0.000 

7.215 

0.9 

0.000 

3.798 

0.000 

4.869 

1.0 

3.^25 

2.^25 

3.971 

2.971 

1.1 

2.578 

1.578 

2.566 

1 .  566 

1.2 

1.6^8 

0.6^8 

1.637 

0.657 

1.5 

1.185 

0.185 

1.095 

0.095 

lA 

0.922 

-0.078 

0.821 

-0.179 

1.5 

0.795 

-0.205 

0.714 

-0.286 

1.6 

0.758 

-0.2^2 

0.709 

-0.291 

1.7 

0.782 

-0.218 

0.767 

-0.255 

1.8 

0.851 

-0.1^9 

0.867 

-0.155 

1.9 

0.951 

-0.0^9 

0.997 

-0.005 

2.0 

1.077 

0.077 

1.148 

0.148 

2.1 

1.167 

0.167 

1.254 

0.254 

2.2 

1.132 

0.132 

1.145 

0.145 

2.3 

1.058 

0.058 

1.025 

0.025 

2.^ 

0.992 

-0.008 

0.941 

-0.059 

2.5 

0.951 

-0.049 

0.906 

-0.094 

2.6 

0.937 

-0.063 

0.912 

-0.088 

2.8 

0.968 

-0.032 

0.986 

-0.014 

3.0 

1.021 

0.021 

1.051 

0.051 

3.2 

1.031 

0.031 

1.050 

0.050 

3.^ 

1.003 

0,003 

0.979 

-0.021 

3.6 

0.985 

-0.015 

0.978 

-0.022 

3.8 

0.986 

-0.01'4- 

1.002 

0.002 
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TABLE  D.7 

THE  COMPUTED  FUNCTIONS  g  AND  H  FOR  THE  HARD  SPHERE 

POTENTIAL 

n*  -  0.9 

r/a  g  H 

^__ 

0.2 

0.5 
0.4 

0.5 
0.6 
0.7 
0.8 
0.9 
1.0 
1.1 
1.2 
1.3 
1.4 
1.5 
1.6 

1.7 
1.8 
1.9 
2.0 
2.1 
2.2 

2.5 
2.4 

2.5 
2.6 
2.8 
3.0 
5.2 
5.^ 
5.6 
5.8 
4.0 


0.000 

5?.S1^ 

0.000 

49.620 

0.000 

42.054 

0.000 

54.970 

0.000 

28.294 

0.000 

22.104 

0.000 

16.495 

0.000 

11.575 

0.000 

7.445 

5-195 

4.195 

2.915 

1.915 

1.550 

0.550 

0.881 

-0.119 

0.626 

-0.574 

0.582 

-0.418 

0.647 

-0.555 

0.774 

-0.226 

0.955 

-0.065 

1.111 

0.111 

1.295 

0.295 

1.556 

0.556 

1.108 

0.108 

0.911 

-0.089 

0.824 

-0.176 

0.829 

-0.171 

0.891 

-0.109 

1.056 

0.056 

1.115 

0.115 

0.992 

-0.008 

0.916 

-0.084 

0.980 

-0.020 

1.052 

0.052 

1.036 

0.056 
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